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ABSTRACT 

The  infinite  Whittaker  summation  and  Shannon's 
sampling  theorem  both  use  the  weighted  sum  of  sine  functions 
(the  "Cardinal"  function)  in  the  interpolation  algorithm. 

When  the  number  of  original  samples  is  approximately  equal 
to  twice  the  product  of  duration  (T)  and  bandwidth  (W) , and 
when  it  is  desired  to  increase  the  number  of  samples  by 
powers  of  2,  the  interpolation  process  can  be  written  as  a 
matrix  equation.  It  is  shown  that  when  the  original 
sample  set  is  periodic,  the  matrix  elements  converge  to 
simple  cosecant  and  cotangent  fxinctions.  ^ 

The  "V'niittaker " matrix  developed  for  either  the  2TW 
transient  or  the  periodic  sample  sets  can  be  manipulated  into 
a real  symmetric  matrix  format.  It  is  shown  that  a unitary 
equivalence  transformation  on  the  periodic  matrix  imple- 
mented via  the  Fast  Fourier  Transform  (FFT)  and  an  orthogonal 

t 

similarity  transformation  on  the  symmetric  matrix  are  really  I 

equivalent  algorithms.  It  is  also  shown  that  by  suitably  | 

modifying  the  transient  Whittaker  matrix,  an  orthogonal  | 


I 

< 

I 

similarity  transformation  is  possible  which  can  significantly 

reduce  the  niamber  of  computations  necessary  to  interpolate  < 

a data  set.  When  "a-priori"  knowledge  about  a data  set  is 

not  available,  but  it  is  desired  to  fit  a specific  curve  to 

the  data,  it  is  shown  that  the  eigenvectors  of  the  modified 

matrix  are  a unique  orthogonal  basis  for  the  space  which 

includes  the  data  set;  linear  combinations  of  the  basis 

vectors  can  then  be  made  using  a technique  called  the  Eigen- 

Filter,  Cross-Correlation  Algorithm  to  indicate  how  well  the 

combination  works.  Numerous  examples  are  given  which  show 

that  this  algorithm  can  provide  a better  interpolation  than 

Fourier  techniques. 

Matrix  norms  and  condition  numbers  are  used  to  bound 
truncation  errors,  computer  rovmd-off  errors,  and  errors  due 
to  the  curve  fitting  algorithm.  It  is  also  shown  that  when 
"noisy"  data  is  interpolated,  the  noise-to-signal  ratio  of 
the  interpolants  can  be  magnified  by  the  interpolating 
matrix  condition  number . 

— ^ An  extensive  computer  program  which  implements  the 
algorithms  is  described.  Numerous  signals  are  processed 
and  the  results  presented  in  plots  and  tabular  form.  The 
work  is  ended  with  an  entire  chapter  suggesting  areas  for 
follow-on  work. 
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INTRODUCTION 

Interpolation  is  a well  established  branch  of  mathe- 
matics. For  the  engineer,  however,  the  tools  to  implement 
interpolation  are  sometimes  lacking;  we  try  to  get  along  with 
simple,  linear  approximations  because  of  the  low  cost  and 
ease  of  implementation.  Because  today's  mass  data  acquisition 
systems  dictate  minimum  sampling  rates,  these  "straight  line" 
approximations  are  often  inadequate.  It  is  our  purpose,  then, 
to  offer  an  alternative.  We  present  algorithms  which  can  be 
implemented  on  big  or  small  computers,  as  well  as  in  hardware; 
but,  as  opposed  to  linear  "approximations,"  our  algorithms 
are  "exact"  when  the  original  phenomenon  and  its  sample  set 
satisfy  certain  requirements. 

Generally,  interpolating  algorithms  are  concerned  with 
fitting  an  "N^^"  degree  polynomial  through  a set  of  equally 
or  unequally  spaced  saimple  points  from  some  continuous 
process.  It  is  well  known  theoretically,  that  N + 1 sample 
points  can  be  connected  by  the  curves  which  plot  all  poly- 
nomials of  degree  N or  greater.  E.  T.  IVhittaker  set  about 
finding  out  if  any  one  of  all  the  fiinctions  that  can  be  made 
to  fit  a data  set  had  any  distinguishing  properties — "a 
function  of  royal  blood  whose  distinguished  properties  set 
it  apart  from  its  bourgeois  brethern."  The  well  known 
cardinal  function  (sum  of  weighted  sine  functions)  resulted 


! 

! 
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from  his  work  and  has  the  property  that  it  contains  no 
frequencies  higher  than  twice  the  sampling  frequency.  In 
more  recent  years,  the  independent  but  equivalent  Shannon 
sampling  theorem  was  developed  and  is  the  basis  for  much  of 
today's  communication  work.  This  theory  also  derives  the 
Siam  of  weighted  sine  ftinctions,  but  as  the  consequence  of 
an  ideal  filter  acting  on  a sequence  of  data  samples  from  a 
band-limited  signal. 

Our  approach  to  interpolation  is  to  formulate  the  sum 
of  weighted  sine  functions  into  a matrix-vector  product. 

We  show  that  the  matrix  with  elements  formed  from  the 
Whittaker  cardinal  function  can  be  written  as  a symmetric 
matrix  with  many  useful  properties.  These  properties  are 
what  we  exploit.  The  rationale  for  this  approach  is  that 
matrices  are  physical  entities;  any  manipulation  of  these 
arrays  is  directly  translatable  to  computer  language  or 
hardware . 

1 . 1 Organization  of  Work 

The  remainder  of  our  work  is  divided  into  several 
chapters  wherein  various  aspects  of  the  interpolation  problem 
are  discussed.  Each  chapter  is  complete  in  that  any 
derivations  begin  and  end  there;  the  Appendix  is  reserved 
for  the  various  computer  programs  rather  than  more  detailed 
derivations  which  this  author  can  better  handle  in  the  main 
text.  Symbology  is  sometimes  complicated,  so  a list  of 


T 


2 


symbols  is  provided  at  the  beginning  of  this  dissertation  to 
aid  the  reader. 

We  begin  in  Chapter  II  by  reviewing  the  mathematical 
process  of  interpolation  in  terms  of  the  digital  filter. 

The  engineering  "band-limited"  function  is  introduced  and 
subsequently,  its  sampled  data  version  is  generated  by  first 
sampling  with  a sequence  of  delta  functions,  and  then  with  a 
more  realistic  sequence  of  rectangular  waveforms.  It  is  this 
sample  data  set  that  we  wish  to  interpolate.  We  show  that  if 
the  samples  are  passed  through  an  appropriate  low-pass 
digital  filter,  the  resulting  continuous  waveform  is  a re- 
constructed version  of  the  original  function.  Practically, 
a more  dense  data  set  rather  than  a continuous  function  is 
the  object  of  the  interpolation  scheme.  Tv/o  finite  inter- 
polation schemes  are  discussed:  Interpolation  by  the  Discrete 
Fourier  Transform,  and  a time  domain  convolution  of  the 
truncated  Whittaker  cardinal  function  and  the  sample  data 
set. 

In  Chapter  III,  we  present  an  historical  review  of 
interpolation  involving  the  weighted  sum  of  sine  functions. 

We  present  the  work  of  E.  T.  Whittaker  [36]  wherein  he 
developed  the  feimous  cardinal  function,  and  we  follow 
through  with  the  work  of  his  son,  J.  M.  Whittaker  [37] , who 
showed  what  types  of  sample  sets  generate  the  cardinal 
function.  The  communications  sampling  theorem  developed 
independently  from  Whittaker  is  discussed  again  but  now  from 


3 


the  point  of  view  of  how  it  relates  to  Whittaker.  We 
present  brief  reviews  of  several  important  papers  by 
Hartley  [14],  Nyquist  [21],  and  Shannon  [30]  which  emphasize 
the  importance  of  how  many  samples  are  needed  to  generate  the 
cardinal  function. 

The  original  effort  in  this  dissertation  begins  with 
Chapter  IV  wherein  we  select  one  of  many  possible  inter- 
polation intervals  and  formulate  the  interpolating  equations 
as  products  of  matrices  and  sample  data  vectors.  The  choice 
of  the  "mid-point"  interpolation  interval  and  square  inter- 
polating matrix  may  seem  arbitrary,  but  we  show  later  that 
once  this  problem  is  mastered,  the  number  of  samples  can  be 
successively  doubled  by  recursively  applying  the  mid-point 
algorithm.  We  also  introduce  the  concepts  of  periodic  and 
transient  matrix  operations;  i.e.,  we  develop  a special 
matrix  to  implement  the  mid-point  algorithm  for  band-limited, 
periodic,  time  domain  signals,  and  develop  a separate  trun- 
cated matrix  for  band  limited,  time  domain  transients.  The 
two  forms  of  the  I’/hittaker  matrix  are  "massaged"  until  they 
are  symmetric  and  the  matrix  elements  can  be  expressed  in  a 
closed  form  trigonometric  expression  for  the  periodic  case, 
or  as  a sequence  of  terms  from  a truncated  infinite  series 
for  the  transient  case. 

In  Chapter  V we  show  that  a matrix  expression  for  a 
periodic  convolution  process  has  very  special  properties  [13] , 
[15] . Specifically,  the  matrix  is  circulant  and  can  be 
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'decomposed  into  the  product  of  a complex  unitary  matrix  with 
a matrix  of  complex  eigenvalues,  followed  by  a product  with 
the  complex  conjugate  transpose  of  the  unitary  matrix.  The 
eigenvalues  are  determined  from  the  Fourier  Transform  of 
one  of  the  rows  of  the  circulant  matrix,  and  the  elements  of 
the  unitary  matrices  are  themselves  samples  from  the  Fourier 
kernel.  We  explore  this  decomposition  for  the  periodic 
Whittaker  matrix  in  cyclic  form;  we  derive  closed  form  ex- 
pressions for  the  matrix  eigenvalues  and  prove  the  surprising 
result  that  the  even  ordered  periodic  Whittaker  matrix  is 
singular.  We  conclude  the  chapter  by  describing  an  inter- 
polation algorithm  which  uses  the  Fast  Fourier  Transform 
(FFT) . 

Both  the  periodic  and  transient  Whittaker  matrices  in 
real  symmetric  form  are  orthogonally  similar  to  a diagonal 
matrix  of  real  eigenvalues.  In  Chapter  VI  we  describe  a 
computer  technique  based  on  the  Francis  QR  [9],  [20]  algorithm 
to  generate  the  orthogonal  matrices  of  eigenvectors  and  the 
diagonal  matrix  of  eigenvalues.  We  then  show  that  this  de- 
composition can  be  viewed  as  a discrete  cross-correlation 
process,  and  if  certain  properties  are  known  to  be  present 
in  the  sample  set,  significant  savings  can  be  achieved  over 
straightforward  implementation  of  the  matrix  products. 

The  purpose  of  Chapter  VII  is  to  present  error  bounds 
for  the  interpolation  algorithms.  Four  types  of  errors  are 
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discussed:  series  truncation  errors,  errors  due  to  noisy 

data,  machine  round  off  errors,  and  errors  caused  by  the 
Eigen-Filter,  Cross-Correlation  algorithm  described  in 
Chapter  VI.  The  benefits  of  casting  the  interpolation 
problem  as  a matrix  process  become  evident  in  this  chapter 
when  we  discuss  the  error  bounds  in  terms  of  matrix  norms 
and  condition  numbers. 

In  Chapter  VIII  we  present  the  results  of  our  work  in 
the  form  of  plots  and  graphs  of  the  outputs  from  the  various 
algorithms.  We  show  that  the  major  algorithms  as  programmed 
in  the  Appendix  do  work  and  are  practical.  Vie  also  summarize 
our  goals  and  findings  as  a conclusion  to  our  work. 

Chapter  IX  is  a special  chapter  wherein  we  outline 
other  problems  associated  with  interpolation.  First,  because 
of  the  mid-point  interpolation  scheme  and  square  interpolating 
matrix,  N - 1 interpolants  are  actually  computed  and  one 
extrapolant  is  produced.  If  the  interpolants  themselves  are 
interpolated,  N - 2 of  the  original  N points  are  returned  and 
two  extrapolants  produced.  Can  this  process  be  continued? 

The  inverse  interpolation  problem  is  also  discussed.  In 
particular,  given  the  interpolants,  how  do  we  compute  the 
original  data  vector?  This  problem  is  not  straightforward 
when  the  periodic  Whittaker  matrix  is  even  ordered;  i.e., 
the  matrix  is  singular.  We  also  outline  a fast  way  for 
computing  the  derivative  of  the  original  scimple  set.  This 


can  be  accomplished  for  the  periodic  case  via  a decomposition 

of  a modified  VJhittaker  matrix  implemented  with  the  FFT.  '' 

Finally,  we  discuss  a recursive  algorithm  which  would  allow 

interpolating  intervals  other  than  the  mid  point  to  be 

approached. 

The  final  parts  of  the  dissertation  are  a compendiiun 
of  computer  programs  in  the  Appendix,  and  a Bibliography. 


CHAPTER  II 


ENGINEERING  APPROACH  TO  INTERPOLATION 

In  Section  1 of  this  chapter,  we  discuss  two  popular 
sampling  functions  used  on  continuous  waveforms:  first,  we 
present  the  idealized  impulse  train  used  for  theoretical 
work  and  then  we  develop  the  realistic  rectangular  pulse 
function.  We  introduce  the  dual  frequency  - time  relation- 
ships of  these  functions,  and  in  Section  2 we  apply  the 
sampling  functions  to  engineering  "band- limited"  signals  to 
produce  the  sample  data  set.  Next,  we  show  the  periodic 
nature  of  the  Fourier  spectrum  of  this  sample  data  set,  and 
then  we  prove  that  when  such  a spectrum  is  filtered  with  an 
ideal  low-pass  filter,  the  original  continuous  time  domain 
signal  is  returned. 

When  a continuous  signal  is  not  needed  from  the  sample 
set,  then  a more  dense  data  set  may  be  the  object  of  an 
interpolation  scheme.  Sections  3 and  4 are  two  approaches 
to  this  problem.  In  Section  3 we  present  the  Discrete 
Fourier  Transform  (DFT)  interpolator — basically,  a frequency 
domain  technique,  and  in  Section  4 we  present  a time  domain 
convolution  approach  to  interpolation. 

The  material  in  this  chapter  is  covered  in  numerous 
textbooks  and  papers.  It  is  presented  here  for  the  sake  of 
completeness  in  discussing  the  interpolation  problem. 
Particulary  useful  references  for  the  first  two  sections 


T" 


are  books  by  Stanley  [31,  Chapter  3]  and  Brigham  [3,  Chapter 
6].  In  particular,  we  follow  Stanley's  lucid  development  of 
the  sampling  functions  in  Section  1 and  use  Brigham's 
pictorial  approach  to  the  sampling  theorem  and  reconstruction 
in  Sections  2 and  3.  Oppenheim's  book  [23,  Chapter  3, 
problem  21]  and  the  papers  by  Schaefer  [27,  pp.  692-702] , 
Urkowitz  [35,  pp.  146-154],  Oetken  [22,  pp.  301-309], 
Crochiere  [6,  pp.  444-456],  and  Rabiner  [26,  pp.  457-464] 
provided  the  motivation  for  Section  3.  Steam's  book  [32] 
provided  an  overall  reference  for  sampling  and  reconstruction 
and  was  particularly  important  because  it  was  the  text  for 
two  of  this  author's  signal  processing  courses  at  the 
University  of  New  Mexico. 


2 . 1 Sampling  Functions 

The  well  known  complex  Fourier  series  representation 
for  a periodic  function  is  expressed  by  Stanley  [31,  p.  38] 


,Tnnt 


(2-1) 


where 


T .,7nnt 

,2  . T 

= T ^ *<•=>« 

~2 


(2-2) 


The  period  T is  the  interval  over  which  the  function  completes 
one  cycle  of  its  periodicity,  and  the  c^^  are  the  complex 
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Fourier  coefficients  of  the  waveform.  The  distinction 
between  x (t)  and  x(t)  is  that  x(t)  is  one  cycle  of  a 
periodic  train  of  cycles,  i.e., 


(t)  = I x(t  - kT) 
k=-“> 


(2-3) 


with 


x(t)  = 


Xp(t) 


_T  < . < T 

2 " 2 


iti  >1 


(2-4) 


Parallel  to  the  Fourier  series  representation  of  a 
signal  is  its  Fourier  transform 


with 


_.2TTft 

X(f)  = |x(t)e  ^ dt 


“ .2TTft 

x(t)  = /X(f)eJ  df 


(2-5) 


(2-6) 


X(f)  is  known  as  the  Fourier  spectrum  of  x(t).  Of 
particular  interest  here,  is  that  when  x(t)  is  a pulse 
(x(t)  =0,  I t 1 > j ) the  Fourier  coefficients  of  the  periodic 
extension  of  x(t),  x (t) , are  simply  given  by 


X(|) 
°m  T 


(2-7) 
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One  sampling  function  widely  used  in  theoretical  work 
is  the  ideal  infinite  impulse  train  consisting  of  unit  delta 
functions  at  intervals  of  extending  to  ±“.  For  such  a 
waveform,  we  can  write  symbolically 


6^(t)  = I 6(t  - kT^)  (2-8) 

P k4-co  s 


Such  a sampling  function,  being  periodic,  has  a Fourier 
series  expansion.  From  equation  2-6  we  write  for  a single 
delta  function 

T 

2 .2TTft 

X(f)  = / 5(t)e"J  dt  = 1 (2-9) 

T 

~2 


then  from  equation  2-7 


c = iilil  = k. 

Ts 


(2-10) 


and  substituting  in  equation  2-1 


«p(t)  * k-  I 


oTTmt 
• “ T 

J S 


(2-11) 


s m=-<» 


The  Fourier  spectrum  of  equation  2-11  itself  is  given  by 


OO  _ 00 


,2TT(f  - ^)t 


1 m''- 

^ I e ^ dt 


(2-12) 


-<*>  s m=-<=° 


which  is  well  known  [3,  p.  22] 
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A(f) 


(2-13) 


= ^ I 6(f  - 
s m=-o°  s 


In  other  words,  the  spectr\am  of  a periodic  impulse  sampling 
train  is  itself  a periodic  impulse  train. 

In  the  real  world,  the  ideal  impulse  train  cannot  be 
generated.  However,  the  rectangular  pulse  train  can  be 
realized.  Proceeding  as  we  did  for  the  impulse  train,  we 
write  (using  symbolic  notation) 


_j;Ut)=  I iQ_iT  - kT^) 


P k=- 
Then,  for  a single  pulse 


” 1 _-;2Trft, 

sinc(f)  = j a (t)  e ^ 


sinirf  a 


TTf  a 


The  Fourier  coefficients  become 


‘^m  T 


. , Tima, 

sin('T=-) 


.Trmg 

' rp  ) 


and  the  Fourier  series  is 


. , Trma , 

sin(-^) 


>Trmt 


— m--»  fILiSS: 

n S m=-«o  I 


.2- 

e^  "s 


s _j  T 


/ n iun  X 

' fp  / 

S 
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(2-14) 


(2-15) 


(2-16) 


(2-17) 


(2-18) 


T 


The  Fourier  spectrum  of  2-18  is  then  found  as 

[m 
T 


sinc(^)  = / ^ I 


sin  ( rp  ) 


m 


'T  ' T ^ ,TTma. 

s -<*>  s m=-“  {--■!  ) 


® e J s dt 


sinc(^) 

s 


1 ” 

= — y 

m L 


„ • m ct\ 

sin  (-;p-) 

s 


T ^ ,iT  m a. 
s m=-<»  (-^) 

s 


6(f  - 


— ) 
T ' 
s 


(2-19) 


(2-20) 


Again,  the  spectriam  of  the  rectangular  pulse  train  is  the 
ideal  periodic  impulse  train  but  amplitude  modulated  by  the 
sinx/x  fxinction. 

This  process  of  postulating  a sampling  function, 
writing  its  Fourier  series  expansion,  and  finding  the 
spectrum,  could  be  carried  out  for  most  any  sampling  wave- 
form. In  any  case,  we  should  arrive  at  the  form  for  the 
Fourier  series  and  spectrum 


.2 


umt 


^ s m=-“  s 


(2-21) 


where 


s If-T- 


T 

2 

/ x(t)e 
2 


-D 


, 27rf  t 


i=s- 


and 


^ I X(^)6(f  - ^) 


L * 

s s 


(2-22) 


(2-23) 
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These  important  properties  of  sampling  functions  (equations  21 
and  23)  are  used  in  the  next  section. 

2 . 2 Generating  the  Sample  Data  Set 

One  of  the  fundamental  properties  of  spectral  analysis 
is  the  duality  of  time  domain  multiplication  and  frequency 
(spectral)  domain  convolution;  i.e., 

00 

h(t)  • g(t)  <=>  H(f)*G(f)  = / H(u)G(f  - u)du  (2-24) 

— oo 

where  "<s=>"  implies  two-way  equivalence  or  duality.  This 
property  allows  us  to  determine  the  effects  of  sampling  a 
continuous  time  domain  signal  g(t)  by  multiplying  by  a 
sampling  function  of  the  form  of -equation  2-21 

00  . 2^-^ 

git)  • X (t)  = |-  I g(t)X(|-)e^  "^s  (2-25) 

^ s m=-®  s 

But,  from  2-23  and  2-24,  this  is  equivalent  to  the  convolution 
00  00  00 

/ G(u)X  (f-u)du=i-  / G(u)  I X(^)5(f-^-u)du  (2-26) 

-oo  m=-“>  "^s  s 

= ^ I X(^)G(f  - |-)  (2-27) 

s m=-“>  ■‘’s  S 

Simply  stated,  the  spectrum  of  the  sampled  waveform  is  the 

periodic  extension  of  the  spectrum  of  the  original  continuous 

waveform,  repeated  at  intervals  of  lm|  =0,  1,  ...  “. 

■‘■s 

These  ideas  are  summarized  in  Figure  2-1  where  we  show  an 
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infinite  cosine  wave  sampled  by  the  ideal  impulse  train. 

The  left  side  of  the  figure  shows  the  time  domain  operations 
The  concept  of  "band-limitedness"  is  naturally  intro- 
duced in  Figure  2-2 (a)  where  we  observe  that  G(f)  is  no 
longer  zero  for  jf | > f„.  The  solid  curve  for  G (f)'  in 

D 3i 

Figure  2-2c  shows  that  because  G(f)  is  not  band“limited, 

G, (f)  is  the  sum  of  the  original  spectrum  G(f)  plus  the 

a 

"tails"  from  all  the  shifted  versions  of  G(f) . This 
phenomenon,  of  course,  is  given  the  name  "aliasing." 

Real  signals  begin  and  end  in  finite  time  and  real 
sampling  schemes  are  of  finite  duration.  These  concepts  are 
introduced  with  the  example  shown  in  Figure  2-3.  Suppose  we 
theoretically  sample  this  g(t)  with  an  ideal  impulse  train 
and  then  convert  to  a "pseudo"  real  world  digital  process 
by  following  the  sampling  with  a rectangular  window  function 
The  resulting  sample  data  set  is  shown  i.  Figure  2-3  (e),  and 
can  be  expressed  as 

g(kT  )-a)(t)  = I [g(t)  •6(t  - kT  ) ]_r;;n_it+^) 
k=-<»  ® 

N-1 

= I g(t) •5^(t  - kT^)  (2-29) 

k=0  P ® 

which  has  the  equivalent  Fourier  amplitude  spectrum  shown 
in  Figure  2-3 (e) 
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|G^(f)*W(f) 

a 


, » _ sin(TTfNT  ) 

I r p /£  ^ \ if s 

1^  L (irfNT  ) 

S s s 


(2-30) 


Any  sampling  function  could  have  been  used  in  place  of  the 
ideal  impulse  train.  The  effects  on  Figure  2-3  would  be 
that  the  amplitude  of  the  repeated,  modified,  spectrum  in 
2-3 (e)  would  vary  as  the  amplitude  of  the  modulating  function 
in  equation  2-27. 

The  intent  of  an  interpolation  scheme  is  to  reconstruct 
the  original  waveform  from  the  sample  set.  Obviously,  if 
we  multiply  the  spectriam  in  Figure  2-1  (c)  by  the  ideal  low- 
pass  filter  response 


F(f)  = 4j47Tl(f) 


(2-31) 


we  get  back  G(f)  in  Figure  2-1 (a).  This  process  is 
equivalent  to  passing  g(kTg)  through  the  ideal  low-pass 
filter  which  returns  g(t),  the  original,  infinite  duration 
time  domain  signal. 

The  problem  is  not  quite  so  simple  in  the  case  of 
Figure  2-2 (c)  where  multiplication  of  the  spectrum  by  the 
ideal  low-pass  filter  response  returns  a corrupted  version 
of  G(f) . The  equivalent  time  domain  process  of  passing  the 
finite  length  data  set  through  the  ideal  filter  returns  a 
corrupted  version  of  g(t) . To  overcome  this  problem,  we 
must  increase  the  sampling  rate  to  conform  with  the  sampling 
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theorem — samples  must  be  taken  at  a rate  at  least  twice 
the  highest  frequency  in  g(t),  or  at  least  a rate  such  that 
aliasing  is  negligible. 

In  the  case  of  Figure  2-3  (e),  multiplication  of  the 

spectrum  by  the  ideal  filter  response  would  have  to  be 

followed  by  a deconvolution  process  to  remove  the  effects 

of  the  windowing  operation.  Alternatively,  we  might  try 

and  choose  windows  which  have  negligible  effects  on  G (f ) . 

sl 

Then  the  ideal  filter  essentially  returns  the  original 
function. 

2 . 3 DFT  Interpolation 

Implied  in  all  the  derivations  thus  far  is  the 
"continuous"  nature  of  time  and  frequency.  But  this  is  not 
quite  satisfactory  to  explain  data  manipulation  on  a 
computer — we  need  a completely  digital  (discrete)  version 
of  the  dual,  frequency-time  relationship  of  discrete  data 
sets.  This  we  provide  by  extending  Figure  2-3 (e). 

First,  sample  the  spectrum  in  Figure  2-3 (e)  (repeated 

»■ 

in  Figure  2-4 (a) ) with  an  infinite  impulse  train  with 
pulses  spaced  l/NT^  apart.  The  equivalent  time  domain 
process  is  convolution  with  another  impulse  train,  and  the 
overall  results  are  two  discrete  periodic  sequences. 

This  process  can  be  written  as  follows: 

N-1 

g(kT^)-a)(t)  = I g(t)6(t  - kT  ) (2-32) 

® k*0  ® 
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® N-1  .27rft 

G (f)*W(f)  = / I g(t)6(t  - kT^)e~^  dt  (2-33) 
-00  k=o  ^ 


N-1  .2TTfkT 

= I g(kT  )e"^  ® (2-34) 

k=0  ® 


Evaluate  at  f = and  we  have  the  periodic  Discrete 

s 

Fourier  Transform  (DFT)  shown  in  Figure  2-4 (c).  Now  we 
write  for  the  fundamental  period 

(2-35) 

N-1  _-2^ 

= I g(kT  )e  ^ , m = 0,  1,  . . . N - 1 

k=0  ® 


We  can  also  write  the  Fourier  series  for  the  periodic  time 
doirain  sample  set  using  equations  2-1  and  2-7 

00  G . 

[g(kTg)  •a)(t)  ]*A(t)  = I ^ ^ (2-36) 

m=-oo  s 

Then,  the  fundamental  period  is 


^k  NT 


N-1^  2lig5i 

I k = 0,  1, 


N-1 


s m=0 


(2-37) 


Equations  2-35  and  2-37  are  a Discrete  Fourier 

Transform  pair  except  for  the  constant  ~ in  equation  2-37. 

•‘■s 

This  can  be  seen  by  substituting  2-35  into  2-37  (without 

the  ^)  . 
s 
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N-1  M-1  -2^^ 

%-k  i ( i v''  I 

P ^ m=0  k=0  ^ 


(2-38) 


N-1  N-1  .2 


7T  (k-p) m 


•«  J. 

= I Z z 9 

r\  1 /\ 


N 


(2-39) 


m=0  k=0 


1 

N m=O^P  " ^P 


(2-40) 


The  interpolation  problem  requires  that  we  increase 
N in  the  above  equations.  To  increase  N by  a factor  of 


2,  we  proceed  as  follows:  Define  a new  sequence 


(g  I = 2k 


h = 


0 = 2k  + 1 


(2-41) 


k=0,  1,  ...N-1 


Then,  equation  2-35  becomes 


2N-1  ^2^ 

F = J f„e^  '^‘^,m=0,  1,  ...,2N-1 
” 1=0  ^ 


(2-42) 


^ N-1 

^m  " Jq  ^2k®  , m - 0,  1 2N  - 1 


N-1  -2^ 

V -IN 


(2-44) 


^m  " ^ ^k® 

k=0  ^ 


, m=0,  1,  ...,  2N-1 


which  has  twice  as  many  spectral  samples  as  does  Now 
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multiply  by  the  ideal  symbolic  filter  function  which 
zeroes  coefficients  from  j to  and  has  gain  2,  and  take 
the  inverse  DFT  of  the  product 


OM  1 II 

_Lyp2j — — I j 2N 


(2-45) 


2 .2^^  2 2N-1  7r£m 

m=|i^+l 


(2-46) 


Let  r = 2N  - m in  the  second  summation;  then, 

^1  2IM!  1 2^^ 

f = i P 2N  ^1  ^p  -j  2N  (2-47) 

„ N m N ^ 2N-r 


N n 

r=2-l 


From  equation  2-44 


p = F * = F 

2N-r  r -r 


(2-48) 


Then, 


N N 

2“^'  27^-  .TTto 

N N , m 


(2-49) 


Finally, 


— y F e- 

N ,N  T . m*" 
m — ( j-1) 


(2-50) 


0,  1,  . . . , 2N  - 1 
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Thus,  midpoint  interpolation  is  achieved  by  computing 
a modified  inverse  DFT  as  prescribed  by  equation  2-50.  The 
zeroing  of  spectral  coefficients  used  to  arrive  at  equation 
2-50  is  equivalent  to  multiplying  the  periodic  spectrum 

by  the  spectrum  of  the  sampled  version  of  an  ideal  lowpass 
filter. 


As  an  excunple,  consider  the  sample  set  generated  by  a 


periodic  sine  wave. 


Then, 


gj^  = sin2-^  , k = 0,  1,  2,  3 


(2-51) 


F„  = I Sin2¥ 

” kio  " 


ty  TTin 

'j  T 


= 0 + e +0-e 


- 


(2-52) 


_.(2m-l)y  __ 

2e~^  sin2^,  m = 0,  1,  . 


..,  7 (2-53) 


Now,  using  equation  2-50 


\ I 


m=-l 


, - , > TT  IT  £.m 

. m„  j-T- 
■'  sin-^  e-"  , 


(2-54) 

— 0,  1,  ...,  7 


_7r  7r£ 

JT  .-r- 


n nZ 


= j {-e^  ^ e“^  ^ + 0 + e e^  ^ } 


(2-55) 


1 -TTZ  -TTZ  _n 

■=•  (je  ^ 4 - je^  4 }5  sin2-^,  £ = 0,  1, 


7 


(2-56) 


25 


which  is  clearly  the  interpolated  version  of  the  original 
sample  set. 


k=-co  il(t  - kT)  (2-59) 

Thus,  g(t)  is  reconstructed  as  a weighted  s\am  of  sine 
functions  with  the  samples  of  g(t)  themselves  serving  as  the 
weights . 


Similarly,  for  Figure  2-3,  we  obtain 


9 

> 


t 


T 


where  G(f)  is  a corrupted  version  of  G(f).  This  is  equivalent 
to 


g (t) 


oo  N-1 

= / ^ g (u) 6 (u-kT  ) 

-«>  k=0 


sin^(t 

•“■s 


U) 

du 

u) 


(2-61) 


sin^(t  - kT  ) 

N-1  ® 

g(t)  = I g(kT  ) — (2-62) 

k=0  ® J-(t  - kT  ) 

■^s 


Clearly  g(t)  can  differ  from  g(t)  depending  on  how  "poorly" 

✓s 

G(f)  approximates  G(f).  As  discussed  before,  this  problem 

can  be  overcome  oy  choosing  u)(t)  such  that  the  convolution 

of  W(f)  and  G (f)  essentially  returns  G (f ) . This  means 

that  G^(f)*W(f)  in  the  fundamental  region  (|f|  < 2^)  is 

■^s 

essentially  the  same  as  G(f).  Then,  when  G,(f)*VJ(f)  is 

a 

multiplied  by  the  ideal  filter  function,  the  resulting  G(f) 
is  also  essentially  the  same  as  G(f) , and  equation  2-62  gives 
a "good"  approximation  to  g(t). 

The  function  expressed  by  equation  2-59  is  formally 
known  as  the  Whittaker  cardinal  function  after  E.  T.  Whittaker 
(1873  to  1956) , the  English  scholar.  Whittaker  arrived  at 
the  interpolating  properties  of  this  function  quite  inde- 
pendently of  the  sampling  theorem  development  used  here. 
Equation  2-59  and  its  time  limited  version  equation  2-62  are 
really  the  beginning  point  for  the  remaining  work  in  this 
dissertation . 
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CHAPTER  III 

AN  HISTORICAL  REVIEW  OF  INTERPOLATING  WITH 
THE  WEIGHTED  SUM  OF  SINC  FUNCTIONS 

In  the  first  section  of  this  chapter  we  briefly  review 
polynomial  interpolation  and  discuss  the  questions 
E.  T.  Whittaker  asked  in  arriving  at  his  cardinal  interpo- 
lating function.  First,  we  show  that  polynomial  interpo- 
lation can  be  a viable  scheme  for  fitting  a continuous  curve 
to  a set  of  data;  then,  we  extend  the  forms  of  polynomials 
to  include  those  generated  by  finite  difference  equations. 
While  studying  the  Newton  difference  form,  Whittaker  was 
bothered  by  the  fact  that  more  than  one  fxanction  could  have 
the  same  finite  difference  table.  He  subsequently  proved 
that  the  cardinal  function  is  the  lowest  frequency  function 
which  passes  through  all  the  data  points  used  to  generate 
the  differences.  We  also  review  the  work  of  Whittaker's  son, 
J.  M.  Whittaker,  who  extended  the  theory  of  the  cardinal 
function  by  specifying  conditions  under  which  it  converges. 

In  Section  2,  we  discuss  the  digital  sampling  theorem  from 
the  point  of  view  of  how  it  relates  to  Whittaker's  theory. 
Particulary  important  is  that  when  the  number  of  samples 
available  for  interpolation  exceeds  the  product  of  signal 
bandwidth  and  its  duration,  the  cardinal  function  and  the 
signal  which  generated  the  samples  are  one  and  the  same. 
Finally,  in  Section  3 we  hint  at  how  sample  sets  might  be 
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manipulated  outside  the  region  of  interest  in  order  to  assure 
convergence  of  the  cardinal  function  inside  the  interval. 

3.1  The  Whittaker  Cardinal  Function 

The  literature  is  replete  with  the  theme  of  polynomial 
interpolation  on  data  sets.  The  ease  with  which  polynomials 
are  generated  and  the  simple  form  of  their  derivatives  and 
integrals  are  no  doubt  fundamental  reasons  for  the  popularity. 
It  is  no  surprise,  then,  that  the  cardinal  function  for 
interpolation  evolved  during  Whittaker's  study  of  polynomial 
interpolation . 

To  provide  a brief  introduction  to  polynomial  inter- 
polation, consider  the  difference  y(x)  - p(x),  where  y(x)  is 
the  given  function  and  p(x)  is  an  degree  polynomial  to 
be  used  to  approximate  y{x).  The  central  idea  in  inter- 
polation is  to  keep  this  difference  small.  Now  suppose  the 
polynomial  takes  on  the  same  values  as  y(x)  at  the  tabular 
points  X = Xg,  Xj^,  ...,  X^.  We  can  anticipate  a result  for 
the  difference  of  the  form  [8,  p.  100] 

N 

y(x)  - p(x)  = R(x  - Xg)  (x  - Xj^)  . . . (x  - x^)  = RII(x) 

(3-1) 

which  is  identically  zero  for  x=x^,  i=0,  1,  ...,N. 

At  emy  nontabular  point  Xj  in  the  interval  Xg  < X^  < X^^  and 
Xj  f Xj^,  we  do  not  expect  this  difference  to  be  zero;  but,  if 
we  define 
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with 


F{Xj)  = y(Xj)  - p(Xj)  - Rn(Xj) 


y(x.)  - p(x.) 


(3-2) 


(3-3) 


n(Xj) 

we  can  force  F(Xj)  to  be  zero.  Now  F(x)  has  at  least  N + 2 
zeroes.  By  Rolle's  theorem  [11,  pp.  61-62] [28,  p.  12],  F'(x) 
is  guaranteed  N + 1 zeroes  between  the  N + 2 zeroes  of  F(x) 
while  F" (x)  is  guaranteed  N zeroes  between  those  of  F'(x). 
Repeatedly  applying  this  theorem  to  equation  3-2,  we  find 
that  (x)  has  at  least  1 zero  in  the  interval  from  x^ 

to  Xj^,  say  at  X = Then,  using  the  fact  that  the  "N  + 1st" 

derivative  of  p(x)  is  zero,  we  can  write 


(N+1) 


(0=0  = - R(N  + D!  (3-4) 


(N+l) 


(N  + 1)  ! 


(3-5) 


Substituting  in  equation  3-1  and  simplifying 


y(Xj)  - p(Xj)  = 


y(N+l)  (5)n(x.) 


(N  + 1)  I 


(3-6) 


Since  Xj  is  any  non-tabular  point  and  since  equation  3-6  is 
true  at  the  tabular  points,  we  replace  Xj  with  x and  write 


y(x)  - p(x)  = I^i")?-''-^' 


(3-7) 


The  behavior  of  equation  3-7  is  difficult  to  analyze. 
However,  it  can  be  shown  [41,  p.  123]  that  if  y(x)  is  an  entire 
function,  (expendible  in  a power  series  which  converges  for  all 
x)  then  the  sequence  of  interpolating  polynomials  Pj^(x)  with 
N = 3,  4,  ...,  defined  on  the  interval  a £ x ^ b,  converges 
uniformly  to  y(x)  on  the  interval.  For  other  types  of  functions 
we  can  say  that  if  some  bound  for  (5)  is  known,  then 

equation  3-7  may  provide  a useful  bound  on  the  error. 

There  are  numerous  forms  for  the  interpolating  poly- 
nomial p(x).  One  widely  used  with  equally  spaced  data  is 
the  Newton  Difference  Formula  [28,  p.  35] 

= Yo  + KAy„  + t (3-3) 

where  the  special  notation  k^^^  is  defined  as 


= K(K  - 1) (K  - 2)...(K  - i + 1) 


(3-9) 


and  A denotes  the  i finite  difference 


Ayg  = Vi  - ^0 


■yg  = A(Ayg)  = Ay^  - Ayg  = 2y^  - y^ 


(3-10) 


A^g  = A^-^i  - A^-^c 


Clearly  equation  3-8  is  true  for  K = 0;  for  K = 1, 

?!  = Yg  + Ayg  = yg  + y^  - yg  = y^ 

for  K = 2 


(3-11) 
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^2  “ ^^0  “ ^0^  '*’  ^2  ~ ^^1  ^0  “ ^2 

(3-12) 

and  inductively  we  can  show  the  values  of  are  cotabular 
with  Yj^.  We  also  note  that  K is  not  restricted  to  integer 
values,  thus  is  defined  at  nontabular  points. 

i\ 

It  was  equation  3-8  that  perplexed  Whittaker  [36,  p.  181] 
He  noted  that  other  functions  have  the  same  difference 
table  as  y(x) 


X 

• 

y 

Ay 

A"y 

Ay 

• 

a-2w 

y-2 

-i. 

• 

« 

a— u) 

y-1 

Ay-1 

y-2 

.2 

a 

a+dJ 

O 1- 

>1  >1 

1 

o 

>1 

<3 

A y_i 

Ayi 

a+2a) 

^2 

• 

1 

A *1 

Ay  ... 


(3-13) 


"...for  we  can  form  a new  function  by  adding  to  y(x)  any 
analytic  function  which  vanishes  for  the  values  a,  a + u), 
a - oj,  ...  of  the  argument,  and  this  new  function  will  have 
precisely  the  same  difference  table  as  y(x)."  He  called  all 
such  analytic  functions  "the  cotabular  set"  and  pointed  out 
that  they  were  all  equdl  at  the  tabular  points  but  in  general 
not  equal  at  nontabular  points.  He  noted  that  "a-priori  " 

f 
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there  is  no  reason  why  p(x)  in  equation  3-8  should  represent 
y(x)  in  preference  to  any  other  function  in  the  cotabular 
set.  Whittaker  then  asked  two  questions:  (1)  "VThich  one 
of  the  functions  of  the  cotabular  set  is  represented  by  the 
expansion?"  (equation  3-8) ; (2)  "Given  any  one  function 
belonging  to  the  cotabular  set,  is  it  possible  to  construct 
...that  function...  which  is  represented  by  the  expansion?" 

In  cuiswer  to  his  questions,  Whittaker  derived 

TT 

“ sin— (x  - a - ntj) 

C(x)  = I f(a  + ntJj)  — ^ (3-14) 

n=-°°  — (x  - a - nw) 

0) 

as  "a  function  which  is  cotabular  with  the  given  function 
y(x),  but  which  has  no  periodic  constituents  of  periods  less 
than  2u)."  He  presented  a lengthy  proof  which  shows  that 
C(x)  is  the  limit  of  p^^  given  in  equation  3-8  as  N goes  to 
“ [36,  pp.  190-192].  The  f (a  + ntJ)  , n = 0,  -1,  ...,  are 
samples  from  any  function  f (x)  in  the  cotabular  set.  The 
fact  that  C(x)  is  generated  from  any  one  of  these  functions 
led  Whittaker  to  call  C(x)  an  invariant  function  - the 
simplest  function  belonging  to  the  set.  He  defined  C(x)  as 
the  cardinal  function. 

Professor  Whittaker's  son,  J.  M.  Whittaker,  made  some 
important  extensions  to  his  father's  work.  Lacking  in  the 
original  work  was  a definitive  statement  of  under  what 
conditions  C(x)  should  reasonably  be  expected  to  converge. 

As  pointed  out  by  J.  M.  Whittaker  [37],  his  father  said  that 


"when  C(x)  is  analyzed  into  periodic  constituents  by  Fourier's 
Integral  Theorem,  all  constituents  of  periods  less  than  2u) 
are  absent; " his  father  then  proceeded  to  produce  an  example 
which  converged  but  which  could  not  be  analyzed  by  Fourier's 
Theorem. 

J.  M.  Whittaker's  results  are  contained  in  his  theorem: 

"If  {fjj)  is  a sequence  of  real  numbers  such  that  the  sum  of 
2 

f^  over  all  N is  convergent,  then  the  cardinal  series  is 
absolutely  convergent  and  its  sum  is  of  the  form 

1 

C(x)  = / { (t)  (t ) cosTixt  + ij;  ( t)  siniTxt  }dt  (3-lb) 

0 

where  4>  and  ip  are  each  square  integrable  on  [0,  1]."  Here, 
the  set  is  the  same  as  implied  in  equation  3-14  except 

we  choose  a = 0 and  o)  = 1. 

Whittaker's  proof  is  fairly  straightforward.  He  notes 
that  due  to  the  Riesz-Fisher  Theorem,  there  are  functions  tj) 
and  ip,  and  a convergent,  square  summable  sequence  such 

that 

1 , 

/ {<{)(t)  - <J)  (t)}^dt  ^ 0 (3-16) 

0 P 

and 

/ {i(;(t)  - ijjp(t)  }^dt  0 (3-17) 


as  p 


“>  when 


= fo 


+ 7 (f  + f 

n -n 


) cosunt 


(3-18) 


n=l 


and 


= 7 (f  + f )sin'rrnt 


n=l 


n -n 


(3-19) 


We  can  show  that  by  multiplying  <J>p(t)  and  4'p(t)  by  cosine 


and  sine  respectively  and  integrating  on  [0,  1]  we  have 


/ ♦ (t)oosiixtdt  = f sisis  + 1 f (f  +£  j ,3ln»(x-n)^SLni.(xi-n), 
- p 0 TTx  2 n -n  Tr(x-n)  Tr(x+n)  ■' 


n=l 


(x+n) 
(3-20) 


and 


'n=l 


(3-21) 


Then 


1 1 
/ ij)  (t)cosTrxtdt  +/  ip  ( t)  siniTxtdt 


sinTr(x  - n) 


n=-p 


"n  ^(x  - n) 


(3-22) 


Equation  3-22  is  the  truncated  version  of  the  cardinal 
function.  Then,  form  the  difference  between  equations 
3-15  and  3-22 
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1 

/ (tjl  (t)  COSTTXt  +Ii/(t)  SinTTXt)  dt 
0 


n=- 


simr(x  - n) 
■n-(x  - n) 


1 1 
/ {4'(t)  - !+)  (t)  }cosiTntdt  - / {'Mt)  - ip  (t)  }simTntdt 

0 P 0 P 


- [/  {(|>(t)  - * (t)}^dt]^/^  + [/  - 4'„(t) 


(3-24) 


which  approaches  zero  uniformly  as  p Thus, 

J.  M.  Whittaker  proved  that  when  a sequence  of  samples  is 
square  summable,  equations  3-14  and  3-15  converge  to  the 
same  function. 

Actually,  much  additional  work  has  been  and  is  being 
done  which  extends  the  cardinal  function  to  ever  increasing 
classes  of  functions.  One  early  work  was  by  J.  M.  Whittaker 
himself  [38]  wherein  he  places  some  additional  restrictions 
on  sample  sets  of  arbitrary  functions  to  gain  convergence  of 
the  cardinal  series.  More  current  work  is  covered  by  McNamee 
[19]  where  the  equivalence  of  the  communication  sampling 
theorem  and  l-Wiittaker ' s theory  is  recognized  and  interpreted 
in  modern  linear  algebra  terminology.  Our  purpose,  though, 
is  not  to  exhaustively  review  the  cardinal  function.  Rather, 
we  have  shown  that  it  is  logical  to  pursue  interpolation  using 
the  VJhittaker  theory. 
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3 . 2 The  Sampling  Theorem  Approach 


Interpolation  with  the  cardinal  function  is  not 
exclusively  in  the  domain  of  the  mathematicians.  Early 
investigators  such  as  Shannon  [29],  derive  the  weighted  sum 
of  sine  functions  completely  independently.  Shannon,  for 
instance,  (page  627)  while  discussing  certain  continuous 
statistical  functions  which  can  be  transmitted  over  a 
communication  system,  writes  an  expansion  for  such  functions 
as 

<=°  ^ simT(2VJx-n) 

S 7r(2Wx-n)  (3-25) 

n=-oo 

He  declares  that  "If  the  function  f(x)  is  limited  to  the  band 
from  0 to  W cycles  per  second,  it  is  completely  determined 
by  giving  its  ordinates  at  a series  of  discrete  points 
spaced  ^ seconds  apart..."  He  notes  that  f(x)  is  represented 
as  a sum  of  orthogonal  functions  with  the  samples  f^^  represent- 
ing the  coordinates  in  an  infinite  dimensional  function  space. 
Furthermore,  "...if  f(x)  is  substantially  limited  to  a period 

T(i.e.,  f = 0,  n>N  and  N = = 2TW)  then  only  2TW 

n X/  zvi 

coordinates  are  non-zero  in  the  function  space.  Thus, 
functions  limited  to  a band  W and  duration  T correspond  to 
points  in  a space  of  2TW  dimensions."  Again  we  have 
Whittalcer's  theorem,  but  with  a different  interpretation; 
while  Whittaker  found  the  lowest  frequency  analytic  function 
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C(x)  cotabular  with  samples  from  any  f(x)  in  the  cotabular 
set  at  the  tabular  points,  Shannon  proves  [30]  that  the 
cardinal  function  and  f(x)  are  the  same  because  band- limited- 
ness and  sample  spacing  criteria  are  met. 

The  real  importance  of  Shannon's  theorem  and  what 
makes  it  so  distinctive  from  Whittaker  is  the  information 
measure  of  "2TW."  Work  earlier  than  Shannon's  by  Lord  Kelvin, 
Hartley,  and  Nyquist  all  sought  or  proposed  similar  quanti- 
fiers for  the  information  content  of  signals.  Modern  communi- 
cation theory  is  based  on  this  important  concept;  i.e., 
we  must  send  and  detect  - as  a minimxim  - these  21^7  coordinates 
of  a signal  space  in  order  to  reconstruct  (interpolate)  the 
original  signal. 

In  his  important  paper  in  1929,  Hartley  [14,  p.  554] 
also  came  to  the  conclusion  "that  the  maximvun  rate  at  which 
information  may  be  transmitted  over  a system  whose  trans- 
mission is  limited  to  frequencies  lying  in  a restricted 
range,  is  proportional  to  the  extent  of  this  frequency  range. 
From  this  it  follows  that  the  total  amount  of  information 
which  may  be  transmitted  over  such  a system  is  proportional 
to  the  product  of  the  frequency  range  which  it  transmits  by 
the  time  during  which  it  is  available  for  the  transmission." 
This  conclusion  was  reached  argumentatively  by  proposing 
various  alternatives  and  rejecting  them  because  they  all 
depended  on  psychological  considerations  as  opposed  to 
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physical  quantities;  e.g.,  the  number  of  symbols  available 
to  a communications  system  should  not  be  used  as  a measure 
of  information  because  if  the  sender  and  receiver  read 
different  languages,  all  messages  could  be  unintelligible. 

About  the  same  time,  1928,  Nyquist  [21,  p.  618]  also 
came  to  a similar  conclusion  while  studying  telegraph  trans- 
mission theory.  In  an  elegant  hueristic  argxament,  Nyquist 
proposes  that  we  consider  an  arbitrary  telegraph  signal  made 
up  of  any  number  and  combination  of  dots  and  dashes  (with  a 
dash  three  times  as  long  as  a dot) ; the  amplitude  of  each 
dot/dash  is  free  to  vary  (the  shape  is  rectangular,  though) 


and  the  message,  whatever  its  length,  is  assumed  to  be  re- 
peated indefinitely  so  that  Fourier  analysis  is  applicable. 
"The  lowest  frequency  component  has  a period  equal  to  the 
period  of  repetition. . .The  next  component  is  double  the 
frequency .. .The  third  component  is  triple  the  frequency,  and 
so  forth.  Certain  components  may  be  lacking. . .while  there 
is  always  a lowest  frequency,  there  generally  is  no  highest. . . 
Next,  Nyquist  supposes  that  we  transmit  such  a signal  and  one 
identical  to  it  except  that  each  element  of  the  new  signal  is 
half  the  duration  of  the  original.  "That  is  to  say,  every- 


thing happens  twice  as  fast  and  the  signals  are  repeated  twice 
as  frequently.  It  will  be  obvious  that  the  analysis  into 
sinusoidal  terms  (by  Fourier  auialysis)  corresponds,  term  for 
term  (with  the  first  signal)  the  difference  being  that  ' 
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next  assumed  that  the  transmission  medium  affected  both 
signals  linearly  (deformed  each  the  same)  and  noted  that 
the  second  signal  would  be  the  exact  counterpart  of  the  first. 
"Generalizing,  it  may  be  concluded  that  for  any  given  deforma- 
tion of  the  received  signal,  the  transmitted  frequency  range 
must  be  increased  in  direct  proportion  to  the  signaling 
speed,  and  the  effect  of  the  system  at  any  corresponding 
frequencies  must  be  the  same.  The  conclusion  is  that  the 
frequency  band  is  directly  proportional  to  the  speed."  In 
other  words,  when  he  doubled  the  speed  he  doubled  the  band- 
width but  cut  the  duration  in  half.  The  information  content 
of  the  signal  was  unchanged,  i.e.,  T • W was  constant  for 
either  signal. 

The  importance  of  the  2TW  concept  to  interpolation  is 
implicit  in  all  our  work  in  this  dissertation.  We  assume 
that  approximately  2TW  samples  are  available;  otherwise  our 
interpolation  schemes  degenerate  to  the  Whittaker  "low 
frequency"  cotabular  algorithm  (aliasing)  which  is  abhorrent 
to  digital  signal  processing.  Perhaps  more  fundamental  in 
our  work  is  that  we  also  view  interpolation  as  a transformation 
or  mapping  of  vectors  in  a linear  vector  space  of  2TW 
dimensions.  The  modern  approach,  then,  is  to  use  matrices 
to  describe  these  transformations  [4,  p.  107],  [10,  p.  32]. 
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3 . 3 A Modified  Sampling  Theorem 


Practical  engineering  phenomena  are  effectively  time 
limited,  i.e.,  they  begin  and  end  in  finite  time  although 
the  precise  instants  may  be  difficult  to  isolate.  A theorem 
by  Paley  and  Wiener  known  as  the  Paley-Wiener  Condition  [24] 
implies  that  such  functions  Ccinnot  be  band-limited  and  thus 
the  concept  of  the  previous  section  is  always  violated:  "A 
necessary  and  sufficient  condition  for  square  integrable 

. (J)  ((d) 

function  A ((d)  1 0 (F((d)  = A((D)e-^  ) to  be  the  Fourier 

spectrum  of  causal  function  (f(t)  **==!>  F(cd))  is  the  convergence 
of  the  integral  " 

«>  I ^■nA(cD)  I 

/ =—  d(D  < “ (3-26) 

-00  1 + (D 

As  used  by  Papoulis  [25,  pp.  219  and  222],  F((d)  = 0,  (Dj^<(d<(D2 
implies  that  A(cd)  is  zero,  and,  therefore,  is  unbounded. 

Thus,  a causal  function  cannot  be  band-limited. 

This  dilemma,  while  troublesome,  is  manageable.  Even 
Shannon  recognized  this:  "...if  f(t)  is  'substantially' 
limited. . .then  only  2TW  coordinates  are  nonzero..."  In 
linear  algebra  terminology,  if  the  coordinates  of  a function 
(signal)  are  essentially  zero  along  all  axes  of  an  infinite 
dimensional  function  space  except  possibly  in  21^^  directions, 
then  these  2TW  coordinates  adequately  describe  the  function. 

A well  known  trick  in  analyzing  nonperiodic  transient 
phenomena  is  to  form  the  periodic  continuation  [8,  p.  417] 
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and  write  a Fourier  series.  It  is  also  well  known  that 
the  manner  in  which  the  function  is  extended  can  have 
serious  consequences  on  the  number  of  terms  in  the  Fourier 
series  expansion;  e.g.,  the  half  cycle  of  a sine  wave  re- 
peated indefinitely  will  have  the  classic  full  wave, 
rectified  sine  wave  expansion  containing  coefficients  out  to 
».  However,  the  simple  artifact  of  repeating  the  half  cycle 
odd  periodically  so  that  a complete  sine  wave  results, 
reduces  the  expansion  to  a single  coefficient.  As  applicable 
to  interpolation  with  the  cardinal  function  the  consequences 
are  obvious  - many  more  samples  are  required  to  interpolate 
the  rectified  sine  wave  than  required  for  the  pure  sine  wave. 

In  his  dissertation,  Campbell  [5]  chooses  samples  from 
finite  duration  functions  which  begin  and  end  with  first 
derivative  discontinuities.  By  forming  the  odd  periodic 
extension  (flipping  about  the  X and  Y axis)  , Caunpbell  derives 
a special  form  of  the  Whittaker  formula  which  capitalizes  on 
the  reduced  bandwidth  of  what  he  calls  the  "regionally 
band-limited  function." 

Our  approach  in  the  next  chapter  does  not  require  such 

restrictions  on  the  sample  set.  We  argue  that  the  simple 

periodic  extension  of  sample  sets  from  practical  engineering 

systems  is  sufficient.  Consider,  for  example,  the  single 

Pa/(/ 

pulse  consisting  of  one  cycle  of  a sine  wave.  By  the  Pauley- 
Wiener  Condition,  its  spectrum  is  infinite.  Therefore,  we 
whould  not  expect  2TW  samples  to  be  exactly  available.  The 
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simple  periodic  extension  of  its  sample  set,  however,  reduces 
the  problem  to  considering  two  samples.  But  our  approach  is 
more  general;  if  we  assume  that  the  digital  sampling  and 
recording  system  is  properly  designed  to  provide  "approximately" 
the  2TW  coordinates  for  all  input  functions  of  interest,  we 
then  have  available  two  different  techniques  for  interpolation; 
first,  the  periodic  extension  of  the  sample  set  can  be  inter- 
polated as  in  Chapter  5;  secondly,  the  transient  sample  set 
itself  can  be  interpolated  as  in  Chapter  6, 
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CHAPTER  IV 


WHITTAKER  INTERPOLATION  AS  A MATRIX  PROCESS 

This  chapter  begins  by  formulating  a matrix  equation 
which  interpolates  one  point  between  every  two  data  points 
in  a finite  length  sample  set.  The  idea  was  struck  upon 
after  reading  Kun-Shan  Lin's  dissertation  [17,  p.  25]  wherein 
he  observed,  in  passing,  some  interesting  properties  of  the 
elements  in  such  a matrix  when  generated  from  the  cardinal 
function.  The  special  problem  of  interpolating  periodic 
sample  sets  is  considered  first,  but  we  then  show  that  when 
the  summation  parameters  in  the  matrix  element  generating 
equations  are  varied,  the  partial  periodic  and  transient 
sample  sets  are  also  interpolated.  When  an  infinitely 
periodic  band-limited  sample  set  is  interpolated,  we  show 
that  the  matrix  elements  can  be  expressed  in  a closed  form 
by  the  cotangent  function.  We  conclude  the  chapter  by 
rearranging  the  matrix  interpolating  equations  into  a 
symmetric  matrix  format.  The  special  properties  of  these 
symmetric  matrices  are  exploited  in  subsequent  chapters. 

4 . 1 The  Whittaker  Matrix  for  Mid-Point  Reconstruction 

Consider  the  problem  of  reconstructing  points  midway 
between  every  two  samples  in  a periodic  scimple  set.  In 
particular,  consider  the  vector  of  N interpolants  generated 
by  a vector  of  M • N original  data  points.  M is  an  odd 
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number  of  periods  (windows) , each  N data  points  in  length. 
In  matrix  notation  we  can  write 


„ M-N 
f = Wx 


(4-1) 


rN  , 


where  f is  the  vector  of  N interpolants  at  intervals  of 
(2i  - l)T/2,  i = 1,  2,  ...N;  is  the  vector  of  M • N 

original  periodic  data  points  at  intervals  of  (j  - 1 + rN)T, 
j = 1,  2,  N,  r = -(M  - l)/2. . .0. . . (M  - l)/2;  W is  the 

N X M • N matrix  of  Whittaker  coefficients 


W = I I 

I I 


(4-2) 


where 


i = 1,  2,  . 


. . , N ; j = 1 , 2 , 


. . , N 

(4-3) 


and 


(J . 


r _ simT[i  - j - rN  + 1/2]  _ ^ (-1 


i-j-rN 


13 


TT[i  - j - rN  + 1/2] 


IT  2(i  - j - rN)+l 
(4-4) 


Figure  4-1  demonstrates  the  data  format  and  partitioning 
for  the  case  of  M = 3 and  a periodic  sine  wave.  Formats  for 
N both  even  and  odd  are  shown. 

Equation  4-1  can  be  simplified  using  equation  4-2  and 
becomes 
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For  M = 

: 3,  N = 

4,  equation  4-5  yields 

f 

~^l-i 

if  1+1 

111 

_i^ 

hif 

1 

^1 

■7-^9 

5 3 11 

I" 

13 

f 

1 1^1 

^1-i 

-i+l+i 

1 

1 

1 

2 

"ll  r 5 

9^1  7 

^1+9 

5 

■3“ 

11 

^2 

1 .1  1 

1 1.1 

1 , 1 

1 

1 

^3 

TT 

13  5 3 

11  3 5 

9 7 

-=r+l+ 

9 

^4 

15 

^1  1 
13  5 3 

1 1^1 

11  3 5 

1 

7 

and 

for  M = 3, 

N = 3, 

■^r 

-fi+i 

-i-i+i 

3 T9 

^2 

_2 

IT 

111 

9 3 3 

i+l-i 

^2 

f3. 

-- — fi+1 
11  5 

111 

9 3 3 

-ifl+i 

y+X+5  ^ 

.^3. 

X, 


X, 


X. 


X, 


(4-6) 


(4-7) 


In  general,  the  summed  periodic  form  of  equation  4-1 


can  be  expressed  as 


where  P is  the  NXN  matrix  with 


and 


P = [Pij]  = I 

r=-(^) 


(4-9) 


i: 


2 

TT 


M-1 

(V> 


/M-lx  Twl 2 


/ I'X  JU  \ 

r=-(— ) 


2{i  - j - rN)  + 1 


(4-10) 


Figure  4-2  is  a plot  of  this  Whittaker  kernel  for  and 

p^^  when  M = 3 and  N = 4,  and  for  p^^^  and  p^^^  when  M = 3 

t h. 

and  N = 3.  Using  equation  4-10,  the  "i  " interpolant  is 
computed  as 


I 


f . 
1 


- N 

- I 

^ j=l 


1 

r=-(!3;i) 


2{i  - j - rN)  + 1 


X . 
3 


(4-11) 


Reversing  the  order  of  summation  and  substitution  of 
j=i+p+l-  rN,  equation  4-11  becomes 


f . 
1 


2 

IT 


' 2 ’ N-i-l+rN 

I I 

_ /M-l>  p=-i+rN 
r“"  \ 2 ' 


i-i-p-l+rN-rN  ^ 


2i-2 (i+p+l-rN) -2rN+l 


i+p+l-rN 


(4-12) 


T 


which  reduces  to 


f . 


1 


= 2 y y (-l)P 

TT  ^ p 2p  + 1 i+p+l-rN 


(4-13) 


til 

Equation  4-11  obtains  from  evaluating  the  "i^  " point 
in  equation  4-8,  while  equation  4-13  is  the  form  for 
evaluating  the  same  point  as  formulated  by  equation  4-1. 

Although  the  two  forms  are  equivalent,  the  number  of  compu- 

t 

tations  is  not  the  same.  For  given  M and  N,  the  bracketed 
terms  in  equation  4-11  can  be  precomputed  and  stored;  thus 

j 

N multiplications  are  implied  to  evaluate  f^.  Equation  4-13,  ’ 

however,  requires  M • N multiplications  for  the  same  i 

computation.  I 

Equations  4-11  and  4-13  can  be  used  when  M is  an  even 
number  of  data  windows  if  the  data  can  be  resegmented  into 
an  odd  number  of  windows  each  containing  an  even  number  of  : 

periods  of  periodic  data  samples.  However,  the  two  important 
cases  are  delineated  by  the  present  form  of  these  equations: 
the  case  where  M = 1,  or  the  transient  case;  and  the  case 

I 

where  M -*•  <»,  or  the  periodic  case.  Incidental  to  these  two  J 

I 

important  extremes  are  the  cases  where  multiple  cycles  of  I 

I 

some  partially  periodic  transient  event  must  be  interpolated.  ' 


4 . 2 Trigonometric  Representation  of  Matrix  Elements 

Equation  4-11  can  be  further  simplified  for  the  periodic 
interpolation  problem.  In  his  dissertation,  Campbell 
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[5,  Appendix  B]  shows  that  the  Whittaker  summation  con- 
verges to  a trigonometric  form  when  the  data  is  periodic 
and  odd  symmetric  about  the  center  of  each  window.  A 
somewhat  different  development,  assuming  no  special  symmetry, 
proceeds  as  follows:  Let  M -*■  and  rewrite  equation  4-10  as 


ID 


= iiil 

ttN 


(-1) 


-rN 


rz=-oo  (i — i + — ) - r 


(4-14) 


If  N is  restricted  to  an  even  number  of  data  points, 
equation  4-14  is  always  positive  and  can  be  written' 


i — n 00 

(-1)  ^ r 1 

^ 


(4-15) 


where 


I 


ii--  3l  + -L 

N 2N 


(4-16) 


The  form  of  the  summation  in  equation  4-15  can  be  observed 
by  writing  a few  terms  around  r = 0 


^ H-r  •••  £+2 


£+2  £+1  £ £-1  £-2 


(4-17) 


Grouping  like  terms. 


y — = 

^ £-r 


(— i—  + - + (— i—  + — — ) + ~ 

'£+2  £-2^  '£+1  £-1^  £ 


(4-13) 


Establishing  a common  denominator  for  each  group  of  terms  and 
simplifying 
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J — = 

^ l-r 


21  . 2Z  . I 

••  ,2_  ^2  ,2  _ ^2  ^ I 


(4-19) 


One  expansion  for  the  cotangent  of  an  argument  [13,  p,  20]  is 


TJl  = — I i 

it£,  it  , „2 


k=i  r - k' 


(4-20) 


Using  equations  4-15,  4-19,  and  4-20,  equation  4-14  can  be 
simplified  to 


(-1)^'^  ±.  , 21  ^ 

N Til  -n  ^ 


o2  2 

I - r 


(4-21) 


(-1) 

Pij  - — 


(4-22) 


t h 

Finally, the  "i  " interpolant  can  be  written  as  the  finite 


length  summation 


N , , . i-j  , 

* jli  -TT  ^ 2n>  Xj 


(4-23) 


For  the  case  where  N is  an  odd  number  of  data  points. 


equation  4-14  alternates  in  sign  and  can  be  rewritten  as 
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j ttN 


I 

k=-00 


(-1) 

l-lk 


-2k 


(-1) 


-2k-l 


i-2k-l 


(4-24) 


_ (-1) 


2ttN 


k=-oo 


S'  1,  1,  1 

2 - I'  2 - - 2 


where  i is  given  by  equation  4-16.  Proceeding  as  before, 
is  simplified  to 


(-1)^"^ 

Pij  = -^2N  [cotTrV2  - cot7T(£  - l)/2) 

and  the  "i^^"  interpolant  for  N odd  is 


(4-25) 


N 


j=i 


^2N^  4N^ 


cot7T(^  + ^ - l)]x. 

(4-26) 


4 . 3 The  Symmetric  Whittaker  Matrix 

We  can  rewrite  equation  4-8  by  reversing  the  sequence 
of  interpolated  points.  Then, 

g^  = Sx^  (4-27) 

N N 

where  g is  f written  in  reverse  order,  and  S is  the 
symmetric  matrix  of  summed  Whittaker  coefficients  obtained 
from  equation  4-10  by  replacing  i with  N + 1 - i 
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N;  j 


If  2|  •••r 


N (4-28) 


s - [S..I  i-  1,  2, 


and 


(_l)N(l-r)  - (i+j)  +1 
2[N(1  - r)  - (i  + j)]  + 3 


(4-29) 


For  the  examples  in  Figure  4-1,  the  interpolation 
formulae  become  for  N even 


.816 

-.058 

-.232 

1.000 

0 

-.652 

^3 

= (.616) 

-.058 

-.232 

1.000 

1.000 

1 

— 

-.759 

^2 

-.232 

1.000 

1.000 

-.232 

0 

.759 

^1 

1.000 

1.000 

-.232 
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-1 

.652 

•• 

_ 

(4-30) 

and  for  N odd. 


■^3“ 

1.049 

-.526 

1.000 

0 

-.889 

^2 

= (.673) 

-.526 

1.000 

1.000 

.866 

= 

0 

^1 

1.000 

1.000 

-.526 

-.866 

.889 

V — 

(4-31) 

The  exact  values  are:  -.707,  -.707,  +.707,  and  +.707  for 
N even;  and  -.886,  0,  and  .866  for  N odd. 

We  can  apply  the  same  procedures  to  equations  4-22 
and  4-25,  yielding  for  the  infinitely  periodic  case  and 
N even 
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(4-32) 


_ __._,2N  - 2i  - 21  +3. 

S.j COt7T( ^ ) 


and  for  N odd. 


‘ i 

I 

I 

I 


i 


i 


(4-33) 


If  we  expand  the  number  of  windows  in  Figure  4-1  to 
the  infinitely  periodic  sine  v/ave,  the  interpolation 
equation  (using  4-32)  becomes 


'"4' 

• 
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1.000 
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= (.604) 
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^1 
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1.000 
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-1 
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(4-34) 

and  using  equation  4-33  becomes 


1.000 

-.500 

1.000 

0 

-.866 

^2 

= (.667) 

-.500 

1.000 

1.000 

.866 

= 

0 

fl_ 

1.000 

1.000 

-.500 

-.866 

.866 

(4-35) 

Figures  4-3  and  4-4  show  the  symmetric  matrix  of 
summed  WhittaJcer  coefficients  (equations  4-29  and  4-33) 
for  N = 9 and  various  choices  for  the  number  of  windows; 


lA 
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Figures  4-5  and  4-6  show  the  matrix  (equations  4-29  and 
4-32)  for  N = 16.  We  note  that  the  matrix  elements  converge 
rapidly  as  the  number  of  windows  increases. 
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4-3  Symmetric  l^ittaker  Matrix 

Nimber  of  Windows  = 1 and  3,  N = 9 
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1*000 

3 0*227 

-0*165  0*1 74  -0*185 

0*  227 

-0*347 

1*000 

1 *000 

-0.347 

4 -0*185 

0*174  -0*185  0*227 

-0  * 347 

1*000 

1*000 

-0*347 

0*227 

5 C * 1 74 

-0*185  0*227  -0*347 

1*000 

1 *0  oc 

-0.347 

0*227 

-0*185 

6 -0*185 

0*227  -0*347  l.^'OO 

1.000 

-0*347 

0*227 

-0*185 

0*174 

7 0*227 

-C.347  1*090  l.OCC 

-<>*  347 

0*227 

-0*185 

0*  1 74 

-0*1 85 

8 -0*347 

1*000  l.COO  -0*347 

0*227 

-0.165 

0*174 

-0*185 

0*227 

9 l.COC 

1*000  -0*347  0*227 

-0*  185 

0*1  74 

-0* 185 

9.227 

-0.347 

4-4  Symmetric  Whittaker  Matrix, 
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Symmetric  VJhittaker  Matrix, 
Number  of  Windows  = 9 and  999,  N 


CHAPTER  V 


CYCLICAL  DECOMPSITION  OF  THE  PERIODIC 
WHITTAKER  MATRIX 

We  begin  this  chapter  by  decomposing  the  periodic 
Whittaker  matrix  into  the  product  of  4 special  matrices 
whose  operations  on  the  data  vector  x can  be  directly  imple- 
mented by  the  Cooley-Tukey  [3],  [34]  factorization  algorithm 
(FFT) . The  first  section  develops  these  matrices  from 
eigenvalue-eigenvector  considerations  and  an  equivalence 
transformation.  In  Section  2 we  rewrite  the  periodic 
Whittaker  matrix  in  a special  cyclic  form.  Finite  difference 
equations  are  then  written  which  express  the  eigenvalue- 
eigenvector  equations  in  discrete  form.  In  Section  3 we 
solve  these  difference  equations  and  obtain  "nice"  analytical 
expressions  for  the  eigenvalues  and  eigenvectors.  A brief 
description  of  the  algorithm  for  implementing  interpolation 
via  the  cyclical  decomposition  is  presented  in  Section  4. 

The  overall  purpose  of  this  chapter  is  to  show  that 

the  matrix  interpolating  equation  4-27  can  be  implemented 

2 

with  fewer  than  the  implied  N operation.  Implicit  in  the 
chapter  is  the  new  interpretation  of  the  workings  of  an 
ideal  digital  filter;  i.e.,  although  we  never  invoke  the 
scimpling  theorem  or  discrete  filter  theory,  we  arrive  at  an 
equivalent  algorithm. 
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5.1 


Equivalence  Transform 

The  eigenvalues  and  eigenvectors  v^  of  equation  4-8 
are  found  from  the  matrix  equation 

Pv^  = A^v^,  0,  i = 0,  1,  N - 1 (5-1) 

The  diagonal  matrix  A of  all  eigenvalues  and  the  corresponding 

matrix  V = [v^  I v^  • • • • I v^  of  all  eigenvectors  also 
•'ll 

solve  equation  5-1 

PV  = VA  (5-2) 


Now  the  symmetric  Whittaker  matrix  equation  4-27  was  developed 
from  P by  premultiplying  with  a permutation  matrix  T which 
interchanges  the  first  and  last  rows,  second  and  next-to- 
last  rows,  and  so  on 

TPV  = TVA  (5-3) 


and  with 


we  have 


S * TP 


SV  « TVA 


(5-4) 


(5-5) 


where  the  elements  of  S are  given  by  equations  4-32  and  4-33. 
If  V is  non-singular,  we  can  rewrite  equation  5-5  as 

S * TVAV"^  (5-6) 

and  the  symmetric  interpolating  equation  4-27  becomes 


62 


g = TVAV'^X 


(5-7) 


Equation  5-6  is  recognizable  as  an  equivalence  trans- 
formation; i.e.,  S and  A are  equivalent  [39,  p.  18].  Our 
.objective  in  this  chapter  is  to  show  that  elements  of  V and 
V ^ are  simply  samples  from  the  Fourier  kernel,  and  that 
the  are  related  to  the  discrete  Fourier  Transform  of 
elements  of  S;  the  entire  equation  5-7  can  be  implemented 
as  a DFT  of  x followed  by  N multiplications  by  the  X^, 
followed  by  another  DFT.  To  accomplish  this,  we  first 
develop  a "cyclical"  representation  for  the  matrix  elements 
of  S,  and  following  the  approach  of  Grey  [13,  p.  17] , we 
solve  the  cyclical  eigenvector  equation  for  A and  V. 

5 . 2 The  Cyclical  Whittaker  Matrix 

Equations  4-32  and  4-33  can  be  rewritten  using  the 
substitution 


£=i+j-2,  0^i+j-2^N-l  (5-8) 


Then , N even 


(-1) 


£+1 


N 


COtTT  ( 


2N 


- 1 - 2£ 
2N 


(5-9) 


N odd 


s„  = 


= (-1) 
2N 


[COtTT  ( 


2N  - 1 - 2£ 
4N 


) - COtTT  ( 


-1  - 2£ 
4N 


) ] 


(5-10) 
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A little  trig  reduces  equations  5-9  and  5-10  for  N even 


s„  = 


= (-1) 
N 


COtiT  ( 


1 + 21. 
2N  ‘ 


and  for  N odd 


s„  = 


= (-1) 
N 


cscir  ( 


1 + 2£, 

2N  ^ 


(5-11) 


(5-12) 


The  corresponding  S matrix  becomes 


where 


®N-2  ®N-1 

®N-1 


®N-3  ®N-2 


= s 


N-X.-1' 


A = 0,  1, 


N - 1 


For  excimple,  for  N = 3 or  4,  we  have 


^ ®2 

=1 

(0 

_f_J 

®1  ®2 

s 

®1  "o 

®2  ®0  ®1 

and 


(5-13) 


(5-14) 


(5-15) 
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t 


k 


* 


y 


®0  ®1  ®2  ®3 

®0  ®1  ^ ^0 

®1  ®2  ^3  ®0 

®1  ®1  ®0  ®0 

®2  ®3  ®1 

®1  ®0  ®0  ®1 

®3  ®0  ®1  ®2 

®0  ^0  ^1  ®1 

L -1 

— — *■ 

The  matrix  equation  5-13  is  of  a cyclic  form  [13], 

[15] , and  [16] . The  system  of  simultaneous  equations  implied 
by  equation  5-5  with  elements  as  shown  in  equation  5-13  can 
be  written  out  for  a general  vector  v and  constant  X. 

®N-2'^N-2'^®N-1^N-1  ^ ^^N-1 
®l'^0^®2''l^®3^2  ■"  •••  ^N-l'^N-2^^0^N-l  = ^^N-2 


®N-3''n-2'''®N-2'^N-1 


= Av, 


(5-17) 


which  is  equivalent  to  the  N difference  equations, 
m=0,  1,  ...,N-1, 


N-1  m-1 

, ^ ®K'^K-m  , ^.®K'^N+K-m  ^'^N-m-l 

k=m  k=0 


(5-18) 


5 . 3 Eigenvalues  and  Eigenvectors  of  the  Cyclic  Form 

One  way  of  finding  the  A and  v for  equation  5-18  is 
to  assume  forms  for  v and  show  that  the  resulting 

expressions  work;  we  assiame  for  the  elements  of  v -f 


65 


(5-19) 


K 


V =0 
K ^ 


Substituting  in  equation  5-18  and  simplifying 


K ^ N K , N -m  -1 

p I Sj^p  + 0 p [ s^p  = Ap  p p 


k=m 


k=0 


(5-20) 


If  we  further  assume 


N 1 
P = 1 


(5-21) 


we  have 


N-1 


* ' I 


K+1 


k=0 


(5-22) 


and  the  general  eigenvector  (normalized)  can  be  written  as 

T 


1 1 2 N-1. 

V = — ( 1,  p,p,...,  p ) 

VN 


(5-23) 


Now  one  function  which  has  the  properties  in  equation  5-19 
and  5-21  is  [1] , [12] 


oils 

p^  = e-^  N 


(5-24) 


Then  we  can  write  for  the  N different  eigenvalues  and  eigen- 
vectors 

^tt(K-H)  I 

A = I s„e"^  ^ , £ = 0,  1,  . N - 1 (5-25) 

k=0  ^ 

^tt£  .Tr£  g (N-1)  tt£ 

£ 1 ,,  -j'^N  -j^N  -j  N 0 r,  1 Ml 

V = -=r-(l,  e , e e ■’  ) , £=0 , 1,  . . . ,N-1 


(5-26) 
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5 

ai 


To  find  the  matrix  of  all  eigenvectors  we  use  equation 
5-26  and  write 


.2^ 

.4^ 

.2-^ 

- 

1 ) IT 

-]  N 

e-7  N 

...  e-^I 

N 

■ 

.8^ 

.4iN- 

- 

Dtt 

-j  N 

e-3  N 

• 

• 

• 

(D 

1 

1 

N 

^(N-Dir  .(N-I)tt 

1 ^ e-^  N • 


(N  - D^tt 


-j“  N 


(5-27) 

which  is  a unitary  matrix;  i,e.,  the  inner  product  of  any 
two  columns  is 


N-i  2^^  2^^ 

.P,  = i I e-^  ^ e^^  N 

K=0 


(5-28) 


,N-1  .2 


y e'^ 

Mi,® 

K=0 


■ff  (p  - q ) K ( 1,  p = a 


0,  p q 


rh«r*fore,  we  can  write  V as  V*,  and 


(5-29) 


1 


1 


1 


1 


V*  = 

VN 


1 e 


1 e 


1 ■5 


. 2- 

. 4- 

2 (N 

- 

1)  TT 

+ 3 N 

+ 1 N 

e ... 

e-^=> 

N 

. 4- 

. 8- 

- 

1)  IT 

+ 3 N 

e+3  N 

e^^ 

N 

, (N-l)  TT 

. (N-l)  TT 

7 (N 

1)^ 

N 

N 

c • • • 

e-^3 

N 

Finally, 


(5-30) 


S = TVAV* 


(5-31) 


where  A is  the  diagonal  matrix  of  eigenvalues  of  P given  by 
equation  5-25. 

We  can  verify  that  equation  5-31  is  true  by  computing 
the  "m,  n^  " element  of  S.  First,  expand  S as 


} 


(5-32) 


with 


^0  0 
V = Tv 


(5-33) 


Then  for  the  m,  n^  element  we  write,  using  equations  5-8 
and  5-25 

N-l  2^  (N  - m)  il  2^  (n  - 1)  £ 

= s „ o=MZe  e ^0  (5-34) 

m,n  m+n-2  ”£=0  ^ 
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fry 


(5-35) 


M-l  2^^  N-1 

I I « I s^e-^  N 

^£=0  K=0  ^ 


tN-1  N-1  .2 

‘•^K=0  ^1=0 


TT (k  - m - n + 2)  £ 


(5-36; 


’m  + n - 2 


Some  further  simplification  for  the  can  be  made. 
When  N is  even,  equation  5-25  can  be  rev/ritten 


.2^  2 ^ -2%^  .2^  N-1  -2%^ 

A,  = e-^  ^ I s N ^ e-^  ^ I s,  e-3  ^ 

^ K=0  ^ , N ^ 

^=2 


(5-38) 


let  p=N  - k - 1 in  the  second  summation,  and  using  equation 


5-14 


.2—  0 


pIT  (p+1)  Z 
J N 


1 

p=5~l 


(5-39) 


ttA  \ |-1  TT(1  t 2k)  £ tt(1  + 2k)il 

- N < r , -j  N +j  N 

= e ^ A 2 Sxf(®  + ® 


(5-40) 


li  |-1  _ 

= e“j  N ^ 2Sj^coSjj-(l  + 2k)  £ 

K=0 


(5-41) 


VJhen  N is  odd,  equation  5-25  becomes 


TT«,  TTkJl  gTT(p+l)  I 

X = e"^  N S y s e’^  N ^ f S ^ + s e’^ 

^ ) kIo  ' P=0  P 

(5-42) 


— 

-jN 


V n TT(1  + 2k)  il 

^ 2s^cos-i— + (-1)  s^ 


(5-43) 


Substituting  equations  5-11  and  5-12  for  we  finally  have 
for  N even 


■ni 

•V  -j  ^ 2,  ,»k  + 2kx ,1  + 2k,  „ ,,, 

Aj^  = e I jj(-l)  cotTr(-2N )cos7r( — ) I (5-44) 

K=0 


Equation  5-44  reveals  that  whenever  ^ = j'  ^ l ~ ^ which  means 
that  even  ordered  periodic  Whittaker  matrix  is  always 
singular.  V7hen  N is  odd,  we  find 

— .,)l+N-l"^ 

^ -j  ^ < V 2,  ,,k  ,l+2k, ,l+2k, „ . (-1)^  2 [ 

An  = e > > 7t(-1)  CSCTT  — ;r7r-)COS7r(-r; ) + r; / 


I jT(-l)  CSCTT  (-^5^)COS7r  (-Tj )Z  + 


(5-45) 


5 . 4 The  FFT  Algorithm 

In  summary,  we  write  equation  5-7  as 


g = (TV)A(V*x) 


(5-46) 


I 
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We  see  that  the  product  V*x  (using  5-30)  is  simply  the 
Discrete  Fourier  Transform  (DFT)  of  the  data  sequence  x. 
Firrthermore,  if  X = V*x,  we  can  write 


and 


F = AX 


g = (TV)F 


(5-47) 


(5-48) 


which  is  simply  an  N point  multiplication  by  the  diagonal 
matrix  of  eigenvalues  with  elements  given  by  equations  5-44 
and  5-45,  followed  by  another  DFT  of  the  product.  The  DFT's 
are  implementable,  of  course,  with  the  Cooley-Tukey  Fast 
Fourier  Transform  algorithm  [34].  Subroutine  "FFTINT"  in 
the  Appendix  implements  this  algorithm. 


CHAPTER  VI 


DECOMPOSITION  OF  THE  REAL  SYMMETRIC 
WHITTAKER  MATRIX 

The  cyclical  decomposition  of  Chapter  V only  applies 
when  the  Whittaker  matrix  S is  periodic.  When  S is  formu- 
lated with  the  transient  generating  equations  in  Chapter  IV, 

S is  still  symmetric  but  has  lost  its  cyclical  properties. 

This  chapter  briefly  reviews  the  orthogonal  similarity  trans- 
formation for  real  symmetric  matrices  as  a prelude  to 
introducing  the  Eigen-Filter,  Cross-Correlation  Algorithm 
which  implements  transient  interpolation.  In  Section  1 we 
prove  that  the  real  symmetric  matrix  S is  orthogonally  similar 
to  a diagonal  matrix,  and  that  similarity  implies  the  diagonal 
matrix  is  the  matrix  of  eigenvalues  of  S.  In  Section  2 we 
discuss  subroutine  SYMEIG  [20]  in  the  Appendix  which  finds 
the  eigenvalues  and  eigenvectors  of  S using  the  Francis  QR 
Algorithm.  In  Section  3 we  show  that  the  similarity  trans- 

I 

formation  on  S can  be  viewed  as  a cross-correlation  process 
wherein  we  first  measure  the  similarity  of  the  eigenvectors 
of  S to  the  data  vector  x before  the  interpolation  process 
is  begun.  By  only  using  significant  correlants,  the  number 
of  operations  necessary  to  implement  g = Sx  can  be  signifi- 
cantly reduced. 

i 
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We  know  that  vectors  in  a linear  vector  space  of 
N = 21^^  dimensions  can  be  expressed  as  linear  combinations 
of  basis  vectors  [10,  p.  75],  We  can  write  arbitrary  vectors 
X and  g with  respect  to  the  basis 


1*2'  ' N 

Q = [q  > q ' • • • ' q ] 
< 1 1 

(6-1) 

as 

1 2 

N 

(6-2) 

X = a^q  + a2q  + . . . 

and 

g = b^q^  + b2q^  + . . . 

(6-3) 

or,  in 

matrix  notation 

X = Qa 

(6-4) 

g = Qb 

(6-5) 

Now,  a 

linear  transformation  S 

which  maps  x to  g 

is  written 

as 

g = Sx 

(6-6) 

Substituting  equations  6-4  and 

6-5 

Qb  = SQa 

(6-7) 

-1 

b = Q -^SQa 

(6-8) 

.( 

! Now  define 
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T = Q"^SQ  (6-9) 

and,  after  solving  for  S,  equation  6-6  becomes 

g = QTQ’^x  (6-10) 

Equation  6-9  is  a similarity  transformation  with  the 
important  property  that  T and  S have  the  same  eigenvalues; 
in  other  words 

det(T  - yD  = det(S  - yl)  (6-11) 


When  S is  a symmetric  matrix,  we  can  write  its 
eigenvalue-eigenvector  equation  as 

Sq^  = Y^q^  (6-12) 

where  i = 1,  2,  . . . , N,  are  the  eigenvalues  of  S and 

the  q^  are  chosen  as  the  associated  eigenvectors.  Using 
inner  product  notation,  we  find 


Yi<q^/  q^>  = <Yiq^,  q^>=<s  q^,  q^>  (6-13) 

or 


By  definition  of  inner  product,  the  denominator  in  equation 
6-14  is  real.  The  numerator  is  real  because  it  equals  its 
own  conjugate.  Therefore,  the  Yj^  are  real  as  the  quotient 


of  two  real  numbers  is  real.  Now,  owing  to  a theorem  by 


T 


Schur  [7,  p.  106],  if  any  matrix  S has  only  real  eigenvalues, 
then  it  is  orthogonally  similar  to  an  upper  triangular  matrix; 
that  is. 


T = Q SQ 


(6-15) 


where  T is  upper  triangular  with  diagonal  elements  t^^  equal 

T T 

to  the  of  S and  with  Q Q = QQ  =1.  This  can  be  shown 
by  hypothesizing  an  eigenvector  for  the  eigenvalue  t^  and 
forming  the  complete  orthonormal  set  Q by  Gram  Schmit  (or 
Householder  as  discussed  later) . Then  we  have 


■qlT- 

■qlT- 

i ^i^^i  = 

2T 

q 

• 

• 

1 * 21 

[Sq-^  1 Sq^  , 
1 1 

1 N 

•••1  Sq"^]  = 

1 

• 

NT 

q 

NT 

q 

1^  2 ^ ^ N 

[t  qtSq^l .* • iSq^l 

■Li  I I 


(6-16) 


t^q  q ' q sq 


I I = |„  .1  (6-17) 

1 


(The  symbol  is  used  to  indicate  that  the  corresponding 
matrix  elements  are  not  germaine  to  the  discussion.) 

We  can  repeat  the  process  on  by  finding  an  ortho- 
normal set  Q2  such  that 
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I 

I 


1 


/ 

t 


After  N - 1 repetitions  we  can  write 


(6-19) 


(6-20) 


T 

But  when  S is  symmetric,  Q SQ  is  symmetric.  Therefore,  T 
must  be  diagonal.  Thus 

S = (6-21) 

where  T is  the  diagonal  matrix  of  all  eigenvalues  of  S. 
Equation  6-21  is  now  an  Orthogonal  Similarity  Transformation. 
Finally,  we  can  write  the  symmetric  Whittalcer  interpolating 
equation  as 
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g = QFq'^x 


(6-22) 


6 , 2 Francis  QR  Algorithm 

Eigenvalues  and  eigenvectors  for  the  symmetric  Whittaker 
matrix  can  in  general  be  computed  from  some  numerical 
technique.  The  method  adopted  in  our  work  is  to  make  this 
computation  using  the  subroutine  SYMEIG  in  the  Appendix. 

This  program  was  developed  at  UNM  by  Dr.  Cleve  Moler  [20, 
Chapter  7]  and  is  based  on  the  Francis  QR  transformation  [9] 

(Q  implies  orthogonal  matrix  and  R is  a right  triangular 
matrix).  Dr.  Holer’s  algorithm  is  the  "Real-Symmetric" 
adaptation  of  the  more  general  technique,  and  consists  of 
two  fiandamental  steps:  (1)  reduction  of  a real- symmetric 
matrix  S to  a tridiagonal  matrix  using  Householder  trans- 
formations; (2)  reduction  of  the  tridiagonal  matrix  to  a 
diagonal  matrix  using  the  iterative  Francis  QR  transformation. 
The  first  step  is  required  because  the  QR  Algorithm  would  be 
too  expensive,  in  terms  of  computer  time,  to  use  on  a 
general  NxN  matrix  S. 

The  first  step  of  the  algorithm  (S  to  tridiagonal)  can 
be  described  in  terms  of  matrix  products;  however,  the 
computer  code  is  quite  different  owing  to  the  fact  that  the 
matrix  operations  simplify  due  to  symmetries  (see  lines  0010 
to  0068  in  SYMEIG,  Appendix) . Given  the  symmetric  matrix 
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Equation  6-27  is  a Householder  reflection  with  the  special 
properties 


P = P 
1 1 


T T 

^ ^1  = ^ 


(6-28) 


(6-29) 


that  is,  Pj^  is  orthogonal.  From  section  6.1  we  know 


S2  = SP^ 


(6-30) 


is  a similarity  transformation  and  that  S2  and  S have  the 
same  eigenvalues. 

The  special  feature  of  the  Householder  transformation 
as  generated  by  equations  6-24  to  6-27  and  as  implemented 
by  equation  6-30  is  the  introduction  of  zeroes  in  the  first 
row  and  column  of  S 


^2  = 


(6-31) 


Now,  form  a P2  which  introduces  zeroes  in  the  first  row  and 


first  coliimn  of  S, 
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The  second  major  step  of  SYMEIG  (tridiagonal  to 


diagonal  via  Francis  QR)  operates  iteratively  on  equation 
6-39.  As  with  the  first  step,  we  describe  the  algorithm 
with  matrix  products  although  the  actual  implementation  is 
different  (see  lines  0069  to  0110  in  SYMEIG) . The  basic 
idea  is  to  let  = Sjj_2  where  Sjj_2  is  as  given  in  equation 
6-39;  then,  factor  A^^  into  the  product  of  an  orthogonal 
matrix  and  a right  triangular  matrix  Next,  reverse 

the  products  and  form  A2  = then  factor  again. 


Continuing,  we  write  the  following  sequence 


Ai 


^2  ~ ^1^1 


A3  - R2Q2  Q3R3 


Ai^  — ^T_  1 Qi  1 — Qi-R. 

k k-1  k-1  k k 


(6-40) 


(6-41) 


(6-42) 


(6-43) 


and  finally. 


*k+l  “ \0k 


(6-44) 


Now, solving  for  Rj^  in  equation  6-4  3 and  substituting  in 


equation  6-44 


^c+1  = Qk^k^k 


(6-45) 


\ 


In  other  words,  equation  6-45  is  again  a similarity  trans- 
formation. The  theory  of  the  QR  Algorithm  is  that  as 
k ->■  oo  the  off  diagonal  elements  of  tend  to  zero;  i.e., 

tends  to  diagonal  form  when  is  tridiagonal.  The 
factorizations  to  are  actually  implemented  by  House- 

houlder  transformations  just  as  we  used  in  going  to  tri- 
diagonal form;  that  is,  we  first  zero  below  diagonal  elements 
in  Aj^  by  N - 2 applications  of  Householder  transformations 

''n-2  "I  \ \ 

then  form 

Qk  = (^^2  ••• 

and  since 


trans  formation 


>■  = ^+1  = 0 so 


(6-50) 


where  F is  the  diagonal  matrix  of  eigenvalues  and  Q is  the 
accumulated  products  of  Householder  transformations.  We 
note  that  Q is  also  the  matrix  of  eigenvectors. 

As  described,  the  QR  algorithm  may  converge  slowly  or 
not  at  all.  However,  it  can  be  shown  [40],  that  by  modifying 
to  Aj^  = A^  - al  where  a is  a root  of  the  lower  2x2  sub- 
matrix of  A^,  and  then  decomposing  A^  as  described  by 
equations  6-40  to  6-45,  that  convergence  is  always  assured. 
Experience  has  shown  that  this  simple  "shift"  generally 
produces  convergence  at  the  average  rate  of  only  one  or  two 
iterations  per  eigenvalue.  This  shift  modification  is 
implemented  in  lines  0083  to  0087  in  subroutine  SYMEIG. 
Figures  6-1  through  6-4  show  the  results  of  applying  SYMEIG 
to  the  Whittaker  matrices  shown  in  Figures  4-1  through  4-4. 


6 . 3 The  Eigenfilter  Cross-Correlation  Algorithm 

If  we  solve  equation  4-29  (element  generating  equation 
for  the  transient  symmetric  Whittaker  matrix)  for 

=1,N'  =i,N-i+2  ^ 2 

these  minor  diagonal  and  minor  subdiagonal  elements  all 

equal  to  2/^.  Then,  the  symmetric  matrix  interpolating 

equation  can  be  modified  to 


g = Sx  = (S  - D) X + Dx  = Sx  + Dx 


(6-51) 
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of  Symmetric 
and  3,  N = 


Now  Dx  is  a simple  linear  interpolating  scheme  and  can  be 
implemented  as  a scaled  stun  of  adjacent  elements  in  x. 

Also,  S remains  symmetric  and  therefore  remains  orthogonally 
similar  to  a diagonal  matrix  of  eigenvalues;  then 

g = Qfo'^x  + Dx  (6-53) 


The  similarity  transformation  can  be  implemented  as 


g 


/N  /N 

QFQ^x  = 


N . 
j=l  => 


x> 


(6-54) 


I 


1 


1 

J 

1 


( 


where  the  inner  product  notation  inside  the  summation 
indicates  that  <q^,  x>  is  a scalar  quantity.  If  we  let 

Pj  = <q^,  x>  (6-55) 


equation  6-54  is  implemented  as 


multiplications  (assuming  Pj^  and  Yj^  are  premultiplied 

together  and  x is  properly  scaled)  and  2N  additions.  3N 

2 

operations  versus  2N  is  very  significant  when  N is  large. 

We  can  interpret  equation  6-60  as  the  sum  of  two 

interpolating  schemes.  The  first  part  is  curve  fitting  in 

a space  of  N dimensions;  i.e.,  when  is  large,  we  say 

the  data  vector  x and  the  eigenvector  q are  very  similar. 

Then,  a part  of  the  interpolated  vector  g consists  of  the 

ttl  ^ ^ Ic 

weighted  version  of  the  "k  " eigenvector,  • The 

second  part  of  g is  given  by  the  linear  interpolating 
scheme  which  consists  of  the  scaled  sum  of  adjacent  data 
points,  Dx. 

Figures  6-5  and  6-6  are  prints  of  the  modified  transient 
Whittaker  matrices  and  the  resulting  eigenvalues  and  eigen- 
vectors when  N = 9 and  16,  respectively.  Figures  6-7  through 
6-10  are  plots  of  the  eigenvectors  when  N = 16.  Subroutine 
"INTERP"  in  the  Appendix  implements  the  Eigen-Filter,  Cross- 
Correlation  interpolating  algorithm. 
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ooooci 


MATRIX* 

00  NMlNOs  0 
♦ 5 

-0*091 
0*111 
-0.1  A3 
0*200 
-0*333 


I / J 

*099 

• 067 
*077 
*C91 

* 1 1 1 

6 -0*143 

7 0.200 

6 -0*333  0*0 

9 0*0  0*0 


0*0 
0*0 

0,0  -0*333 

-0*333  0*200 


6 

0*  1 It  -0*143 
•0*143  0*200 

0*200  -0*333 
0*333  0*0 

0*0  0*0 
0*0  -0*333 

-0*333  0*200 

0*200  -0*143 
0*143  0*111 


7 

8 

0*200 

-0*333 

0*333 

0*0 

0*0 

0*0 

0*0 

-0*333 

0*333 

0*200 

0*200 

-0*143 

C*143 

C*  1 1 1 

0*1  1 1 

-0*091 

0*091 

0*077 

PRINT  CF  MODIFIED 
MULTIPLIER*  C .6376 
2 3 

C.C67  0*377 

0*077  -0*091 
C.C9I  0*111 
0*  111  -0*143 
0*143  0*200 

0*200  -0*333 
C*333  C*0 


9 

0*0 

0*0 

-0*333 

0*200 

-0*143 

0*111 

-0*091 

0*077 

-0*067 


PRINT  OF  eigenvalues* 


MULTIPLIERS  C*982E 

00  NminDs  0 

I/J  1 

2 3 

4 

5 

6 

7 

8 

9 

1 )*9S6 

-l.OOC  0*503 

0*324 

0*129 

-0*2  I 8 

-0*466 

-0.411 

-0*006 

PRINT 

CF  EIGENVECTCRS* 

MULTIPLIERS  0*6C6E 

00  NMINDs  0 

I/J  1 

2 3 

4 

5 

6 

7 

8 

9 

1 0*190 

-C.S09  0*913 

0*  1 77 

-1 *000 

0*598 

0*391 

0*227 

0*943 

2 -0.4C1 

0*462  -0*623 

0*817 

-0*212 

0*589 

0*130 

0*368 

9*867 

3 0.770 

-0.213  0*387 

0*589 

0*708 

0*354 

-0*847 

0*499 

0*044 

4 -0*887 

-0*179  9*384 

e*490 

0*  563 

-0*235 

0*522 

0*726 

-0.589 

9 0*746 

0.S64  -0*382 

0*234 

-0*589 

-0*615 

0*202 

^*  774 

-0*534 

6 -0*389 

-0*812  -0*214 

-0*279 

-0*403 

-0.69C 

-0*654 

0*644 

0*599 

7 -0*058 

C*838  0*654 

-0*  891 

0*237 

0*070 

0*0  56 

0*646 

0.563 

8 C.377 

-0.643  -0*726 

-0*637 

0*374 

9*65  1 

0*624 

0*928 

-0*074 

9 -0*946 

0.285  -0*244 

-0*337 

-0*397 

0*745 

-0*839 

0*232 

-0*338 

6-5  Modified  Transient  Matrix, 

Number  of  Windows  = 0,  N = 9 
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Modified  Transient  Matrix 
Number  of  Windows  = 0,  N ■ 
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Eigenvectors  (5-8)  of  Modified  Matrix 
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CHAPTER  VII 
ERRORS 


The  relative  ease  in  deriving  error  bounds  for  matrix 
equations  was  one  of  our  motivating  factors  in  developing 
the  Whittaker  matrix  process.  We  begin  this  chapter  by 
briefly  reviewing  vector  and  matrix  norms  and  condition 
numbers,  and  then  use  these  concepts  in  deriving  various 
error  bounds  in  Sections  2 through  5.  In  section  2 we 
derive  two  expressions  which  describe  the  effects  of  the 
truncated  Whittaker  matrix  on  the  interpolated  data.  First, 
we  bound  the  error  caused  by  the  truncation  and  then  present 
a formula  which  predicts  the  sensitivity  of  the  interpolants 
to  data  outside  the  interpolating  interval.  In  Section  3 
we  develop  an  expression  for  the  noise-to-signal  ratio  of 
the  interpolated  data.  We  show  simply  that  this  ratio  is 
the  magnified  noise-to-signal  ratio  of  the  original  data. 

The  magnification  factor  is  shown  to  be  the  condition  number 
of  the  Whittaker  matrix.  We  discuss  the  effects  of  the 
computer  hardware  on  the  interpolants  in  Section  4 . By 
assuming  the  worst  possible  computer  roundoff  errors,  we 
show  that  interpolation  with  the  Whittaker  matrices  is 
relatively  immune  to  roundoff  problems.  In  the  final 
section,  we  derive  error  expressions  for  the  Eigen-Filter, 
Cross-Correlation  Algorithm  described  in  Chapter  VI. 
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7 . 1 Norms  and  Condition  Numbers 


The  norm  of  a vector,  signified  by  ( | • | | * is  a 
functional  from  vectors  to  non-negative  reals  which 
satisfies  the  following  properties  [33,  p.  163] 


lx  >0 


Mx|  I = 0,  X = 0 


II  Axl!  = I Al  l|x|| 

X + y I I - 1 |x| I + I |y  1 


(7-1) 

(7-2) 

(7-3) 

(7-4) 


As  used  in  this  chapter,  we  define  the  Holder,  or  "p",  norms 
as  [33,  p.  1661 

N 


Then, 


|x| Ip  = ( I |x  |P)^/P,  1 < p < 

P i=l  ^ 


N 


|x|  Ij^  = I |x^ 

^ i=l  ^ 


(7-5) 


(7-6) 


N 1/2 


|x| I2  = ( I 1x^1  ) 

^ i=l  ^ 


(7-7) 


M M max  I 

I lx|  L = i l*i 


(7-8) 


Given  any  vector  norm,  the  subordinate  matrix  norm  is 
defined  as  [39,  p.  56] 
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r'-.x. 


I 

I 


MAX  I ISx I I 
llx||=l 


(7-9) 


I 

In  addition  to  the  properties  in  equations  7-1  through  7-4 
(replace  x with  S) , matrix  norms  satisfy  the  following 


IISxIl  < llsll  11x11 


(7-10) 


which  follows  from  equation  7-9. 

Some  specific  norms  have  additional  useful  properties. 

T T 

If  S is  symmetric  and  Q is  orthogonal  (i.e.,QQ=QQ  =1) 
we  have  [33,  p.  308] 

llsllj  =||QrQ’'|l2  = Ilr|l2  - lY„axl 

where  the  y ^ are  eigenvalues  of  S,  and  T is  the  diagonal 
matrix  of  all  eigenvalues.  We  also  have  the  "Frobenius" 
norm  [33,  p.  173]  which  can  be  defined  for  a real  symmetric 
matrix  as 

5 1/2  1/2 

llsllp  - sjj)  = '|‘kk' 

where 


T = S^’s  = (QrQ'^)'^(QrQ'^)  = (7-13) 

Then 


and  we  have  the  useful  results 


= Iwl  - I Is  I Ip  = ^ yJT  li,  I 

1 

(7-15) 

The  condition  number  of  a matrix  is  defined  [33,  p.  190] 
as  the  product  of  the  norms  of  the  matrix  and  its  inverse 

= i 1^1  Ip  I IS'^I  Ip  (7-16) 

When  S is  symmetric  and  p = 2,  we  have  using  equation  7-11 

k2(s)  = (7-17) 

^min 

We  note  that  large  implies  "*■  0*  In  other  words, 

"poor"  or  a large  condition  is  predicated  on  "closeness"  to 

2 

singularity.  Also,  when  S = I,  we  have  ^2(3)  = 1.  Thus, 
involutory  or  orthogonal  matrices  are  "perfectly"  conditioned 
with  respect  to  the  2 norm. 

In  the  remainder  of  this  chapter  we  drop  the  subscript 
p from  our  derivations.  Unless  otherwise  indicated,  we  will 
always  restrict  our  discussion  to  forms  involving  the  "2" 
norm. 

7 . 2 Truncation  Errors 

Truncation  errors  result  when  we  truncate  the  infinite 
Whittaker  summation  and  formulate  the  transient  Whittaker 


! 


I 

I 


matrix  interpolating  equation 


g = Sx 


(7-18) 


Then 


and 


X = S ^g 


|xl 


< I I ^ 1 


Mgll 


(7-19) 


(7-20) 


From  equation  4-13  write 


fN  = 2 (-l)P^ 

^i  IT  ^ . 2p  + l*i+p+l 
p=-i 


(7-21) 


where  the  superscript  "N"  indicates  that  we  have  truncated 
the  Whittaker  summation  to  N terms.  Now,  suppose  we  inter- 
polate with  an  N + L point  scheme  and  formulate  the 
difference 


If^^^  - f^l  = i 

' 1 1 ' TT 


N-i-l+L 

Jn-I  2p  l^i+p+l 


(-l)P. 


Define 


AC  T • I jN+L  -N I 

Af . = Lim  f . - f • 

1 T ^ 1 1 


Then,  using  the  Cauchy  Schwartz  inequality 
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(7-22) 


(7-23) 
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< 2 

1/2 

00 

y X ^ 

“ ir 

p-Vi  2P  * l’ 

p=Li  p 

1/2 


(7-24) 


From  section  3.1  we  assume  that  a data  set  {x^,  i = 1,  2,  ...» 
“}  must  be  square  summable  in  order  to  use  Whittaker  inter- 
polation. Thus,  we  assume  the  norm  of  the  data  outside  our 
interpolating  interval  is  equal  to  some  constant  a times  the 
norm  of  the  data  vector  x;  i.e.. 


I ^ ' 

p=N+l  P 


‘ 1/2 

< 

V 2 

- a 

Li=i 

1/2 


= a X 


(7-25) 


Using  this  inequality  and  equation  7-20,  we  have 

1/2 


Af . 5 

1 IT 


( 1 ) ^ 


J«-i  '2P  ^ 1 


,-l 


(7-26) 


Recognizing  that  Af^  = ^'Jjj-i+l'  state  a bound  for 

the  normalized  truncation  error  as 


'N-i+1 


Ag„  ...  <2 
^N-i+1  - 


|5|  1 


I 

p=N-i 


( — ^ ^ 
^2p  + 1^ 


1/2 


or 


< 2 
®i  - 


“ 1-2  - 

1/2 

1 

P=0  P=0  ^ j- 

Y . 
min 

is'^ll 


(7-27) 


(7-28) 


^ ^a[1.2337  .. 


i-2  - 1/2  , 

Iq^2p  + 1 ^7” 

p=0  ^ 'min 


(7-29) 


Easily  computed  for  N = 16,  we  have 


e^^  f |a[1.2337  - 0.0]^^^ 


.385 


= 1.837a 


(7-30) 


*16  - I*I1-2337  - 


1.2170]^'^^-^  = .213a 


(7-31) 


If,  for  instance,  we  assume  a is  .1,  e^  1 .184. 

Equation  7-27  can  be  put  in  a more  convenient  form  if 
we  approximate  the  summation  with  the  integral 


00  00 

L <2rVT>^-^/  . 


p=2N-i  2P 


N-i  (2x  + 1) 


■2°^  ~ 4N  - 4i  +2 


(7-32) 


Then,  the  error  expression  becomes 


^ It-  I — I 

1 IT  Ml  -2^  'y  . ' 


(7-33) 


For  N = 16,  ^ .210a  which  agrees  with  equation  7-31. 

Another  factor  relating  to  truncation  error  is  a 
sensitivity  factor.  First,  we  write  equation  7-18  with 
superscripts  indicating  the  number  of  data  points  used  in 
the  interpolation.  Thus 


N -N  N 
g = S X 


(7-34) 
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If  we  add  one  more  data  point  to  the  interpolation  scheme, 
we  can  write 


N+1 


-N+l 


*N+1 


(N+1)T 


'N+1 


.N+1 

’n+1 


N 


^+1 


(7-35) 


N ^N+1 

Subtracting  g from  g and  multiplying  the  absolute 
difference  by  Ijx^lj  we  can  write 

;N+1  _ N,  I I N,  , < , -N+1  , , ,„-l 


g"  ■ - g 1 1 |x"  I I |s-  I |s 


For  the  i interpolant,  define 


6.  = 

1 


Igr'j  g"' 

-N 


N 


(7-36) 


(7-37) 


then 


r < I N+1 I 'n+1 ' 1 


(7-38) 


If  we  use  equation  4-29  for  the  transient  case 


6.  ^ 2 

1 IT 


Then,  for  i = N 


. < 2 

°N  TT  13  - 2N 


1 1 '^N+l* 

1 

3 - 2i ' II  N 1 1 

1 |x  1 1 

1 min  1 

1 1 1^+1  ^ 

1 

Y, 


min 


(7-39) 


(7-40) 


and  for  N = 16,  we  can  now  state  the  normalized  sensitivity  of 


the  first  of  16  interpolants  to  the  change  caused  by  adding 


one  more  data  point 


6 ^11.  ^+1 
16  n 29  |,  N, 


(7-41) 


^16- 


I 1.  N, 


(7-42) 


Loosely  stated,  if  a 17  data  point  is  about  the  size  of 

the  norm  of  the  first  16  points  used  in  the  interpolation 

(a  rather  horrible  situation)  then  we  should  expect  up  to 

a 6 percent  change  in  the  interpolant  if  we  go  to  a 17 

point  scheme.  However,  if  we  assume  is  only  about  .1 

times  the  norm  of  the  first  16  points  (still  a bad  situation) 

then  the  interpolant  is  relatively  "insensitive"  to  the 

17^^  point;  i.e.,  5 - .006. 

Xd 

7 • 3 Errors  Due  to  Noisy  Data 

Suppose  the  data  vector  in  equation  7-20  is  corrupted 
with  noise.  We  say 


(g  + Ag)  = S (x  + Ax) 


(7-43) 


where  Ax  is  the  vector  of  instantaneous  noise  values  and  Ag 
is  the  resultant  change  in  the  vector  of  interpolants . Now 
we  can  write 


Ag  = SAx 


(7-44) 
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lAgll  1 


lAxll  = llAxI 


(7-45) 


llxll  1 Ms- 


(7-46) 


we  have 


|Ag||  Mx||  < llsll  Ms' 


iAx| I I |g| 


= k(s)  I |Ax|  I I |g| 


(7-47) 


or  finally. 


il.Agl  I ^ k(s)  IjAxI 
M^ll  ^ ll"M 


1^1 

'min 


(7-48) 


where  k(s)  is  given  by  equation  7-16. 

The  ratio  | |Ax| |/| |x| | is  the  ratio  of  the  square  root 
of  noise  power  to  square  root  of  signal  power,  i.e.,  roughly 


the  reciprocal  signal-to-noise  ratio,  (S/N) 


Thus , 


equation  7-48  implies  that  Whittaker  interpolation  can 
magnify  the  N/S  ratio  by  k(s).  For  the  periodic  symmetric 
nonsingular  t^hittaker  matrix  with  N = 9,  k(s)  = 1,  thereby 
preserving  N/S.  For  the  16  point  transient  interpolation, 
k(s)  = 2.6,  which  means  that  the  interpolants  have  gained 
less  than  8.3db  in  noise.  The  value  for  k(s)  for  any 


specific  transient  problem  can  be  computed  by  using  sub- 
routine SYMEIG  in  the  Appendix. 
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7 . 4 Machine  Roundoff  Errors 

We  know,  in  general,  that  the  elements  of  S and  x 
cannot  be  exactly  represented  internally  in  a computer.  We 
can  express  the  resulting  errors  in  g as  follows 

(g  + Ag)  = (S  + AS) (x  + Ax)  (7-49) 

where  AS  and  Ax  are  now  the  machine  errors  caused  by  round- 
off. Then 


(g  + Ag)  = Sx  + [SAx  + AS(x  + Ax)]  (7-50) 


or. 


Ag  = SAx  + AS(x  + 

Ax) 

(7-51) 

Proceeding  as  before,  since 

lUII  i lls-^ll  ||g| 

1 , we  write 

1 |Ag|  1 1 |x| 1 1 1 !sax  + 

4S(x  + Ax)  1 1 1 |S''-|  1 

llg|| 

(7-52) 

and  using  equations  7-4  and 

7-10, 

UMli  < k(s)  + 

itrfii  ' 'liivli 

IIASI.I  ll-xM  ^ hall 

iiCii  iivii  iiCll 

I ilAxii, 

1 1 V i t J 

(7-53) 

On  the  IBM  360/67  using  single  precision  arithmetic, 
roundoff  can  be  expressed  as 
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(7-54) 


= g.(l  + 6^),  |6^|  < 16 

If  we  assume  the  worst  possible  case  for  \ |Ax1 [ and  [ [ASj  j 
we  can  write 

llAxjl  = ||l6”^x||=  16"^  ||x||  (7-55) 

Also, 

I lAS|  I < 1 |AS|  Ij.  = 16’^1  |S|  Ip  (7-56) 

Using  equation  7-15, 

llASlIp  ^ 16'^(lYi)^^^  ^ (7-57) 

But  |y  , I is  simply  l|s||;  therefore, 

lUaX ' 

1 |AS| I < 16"^/n| 1S| 1 (7-58) 

Substituting  in  equation  7-53 

j j I < k(s)  [16"^  + /n  16"^  + VT  16"^°]  (7-59) 

For  moderate  N,  then,  we  can  reasonably  expect 

I I ' ^ k(s)  16”^  = 16“^  (7-60) 

M ’ I I '^min 

For  the  transient  Whittaker  matrix  and  N = 16,  equation  7-60 

implies  that  the  norm  of  the  normalized  error  in  the  inter- 

-5 

polants  due  to  roundoff  is  less  than  10 
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7 . 5 Approximation  Errors 

In  the  Eigen-Filter,  Cross-Correlation  Algorithm,  we 
produce  an  approximation  g*  to  g.  If  we  solve  the  following 
equation  for  x*. 


we  find  that  x - x*  is  in  general  non-zero.  Define  a vector 
of  residuals 


r = X - X* 


(7-62) 


Then, 


r = X - S 


-1 


(7-63) 


and  multiplying  both  sides  by  S,  we  have 

Sr  = Sx  - g*  = g - g*  (7-64) 

But  g - g*  is  the  error  due  to  the  approximation.  Then, 
define  an  error  vector 


e = Sr 


(7-65) 


Proceeding  as  in  previous  sections 


lell  < IISl 


(7-66) 


Using  I |x 


* I I < 


.-1 


|g* I I , we  express 


c*ll  < Ms 


,-l 


IrM  Mg^ 


(7-67) 
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or 


This  "nice"  expression  allows  the  norm  of  the  normalized 
error  in  g*  to  be  bounded  without  really  knowing  g;  i.e., 
compute  g*  as  described  in  Chapter  VI,  use  equations  7-61 
and  7-62  to  find  x*  and  r;  and,  finally,  use  equation  7-68 
to  bound  the  errors  in  g* . For  the  transient  Whittaker 
matrix  and  N = 16,  the  norm  of  normalized  errors  due  to 
approximation  is  less  than  2.6  times  the  norm  of  normalized 
residuals  in  x. 
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CHAPTER  VIII 


RESULTS 

The  Appendix  is  the  listing  of  a lengthy  computer 
program  which  implements  the  major  algorithms  of  this 
dissertation.  The  program  evolved  over  a year's  work  and 
is  somewhat  elaborate  as  to  documentation.  We  believe  that 
reading  Section  1 of  this  chapter,  wherein  we  briefly  describe 
the  program,  and  then  reading  the  listing  comments  on  input 
and  output  parameters,  will  enable  one  to  run  the  program  if 
so  desired.  Section  2 of  this  chapter  is  basically  a 
compendiiam  of  results  from  the  program  and  a discussion  of 
their  significance.  Ten  rather  difficult  to  interpolate 
functions  were  used  to  produce  the  various  figures  in  this 
section.  The  final  Section  3 is  a summary  of  what  we  set  out 
to  do  in  the  dissertation  and  our  own  assessment  of  what  we 
actually  accomplished. 

8 . 1 Description  of  Computer  Program 

The  techniques  discussed  throughout  the  dissertation 
were  programmed  as  docummented  in  the  Appendix.  Basically, 
we  progreimmed  a three  step  effort:  on  the  first  pass,  the 
transient  symmetric  Whittaker  matrix  is  formed  in  subroutine 
MATRIX,  and  subroutine  SYMEIG  is  used  to  find  the  eigen- 
values and  eigenvectors.  Subroutine  SIGNAL  produces  N data 
points  from  some  band-limited  function  and  subroutine  INTERP 
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produces  a vector  of  interpolants  by  implementing  the 
orthogonal  similarity  transformation  of  equation  6-58.  The 
routine  also  produces  a vector  of  correlation  coefficients 
per  equation  6-59.  Various  prints  and  plots  are  produced 
by  subroutines  PRINT  and  PLOT  to  document  the  first  pass. 

Before  the  second  step,  one  analyzes  the  results  of 
the  first  pass  and  selects  eigenvectors  that  "strongly" 
correlated  with  the  original  data  vector.  We  considered  a 
correlation  coefficient  of  .2  or  greater  as  being  significant. 
These  eigenvector  numbers  are  input  for  the  second  step. 

This  run  basically  repeats  the  first  step  except  that  only 
the  selected  eigenvectors  are  used  in  INTERP  to  produce  the 
vector  of  interpolants.  Again,  various  prints  and  plots  are 
produced  to  document  the  second  pass . 

The  third  step  is  for  the  periodic  symmetric  V'Thittaker 
matrix.  Subroutine  MATRIX  produces  this  cyclical  form  of 
the  Whittaker  matrix  and  the  other  subroutines  follow  as  in 
the  transient  case:  subroutine  INTERP  interpolates  using 
the  orthogonal  similarity  transformation,  correlation  co- 
efficients are  produced,  and  prints  and  plots  are  produced. 

In  addition,  subroutine  FFTINT  is  invoked  to  implement  the 
equivalence  transformation  discussed  in  Chapter  V.  Both  the 
INTERP  and  FFTINT  routines  are  used  to  document  that  they 
produce  the  same  results  (except  for  roundoff)  although 
being  completely  different  algorithms. 
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In  the  final  runs  for  this  dissertation,  all  three  steps 
were  run  together  since  prior  experimentation  determined 
which  eigenvectors  were  to  be  used  for  the  second  step. 

The  compile  and  run  times  are  documented  in  Figure  8-11. 

In  addition,  the  program  prints  a total  of  2200  lines  when 
the  program  listing  is  requested,  or  about  1300  lines  for 
all  three  steps  when  only  results  are  required.  Total  cards 
in  the  program  Fortran  deck  are  about  850,  including  comments. 


8.2  Comparison  of  Relative  Errors 

J 

Figures  8-1  to  8-10  depict  the  results  of  running 
the  program  in  the  Appendix  on  10  selected  signals.  (The  \ 

signals  are  shown  as  solid  lines.)  The  "a"  part  of  the  ' 

Figures  are  plots  of  the  interpolants  using  the  transient 

I 

technique  of  Chapter  VI,  and  also  show  the  correlation  co- 
efficients from  equation  6-59.  The  "b"  parts  show  plots  ’ 

of  the  interpolants  using  selected  eigenvectors  and  plots 
of  the  interpolants  from  the  FFT  algorithm  of  Chapter  V. 

The  final  "c"  parts  of  the  Figures  are  tabular  summaries  of  ^ 

the  data  used  to  generate  the  plots.  The  upper  left  table  ^ 

is  for  the  full  transient  interpolation  scheme  and  the  next 

I 

table  is  for  the  selected  eigenvector  approach.  The  lower 

two  tables  are  for  the  symmetric  periodic  matrix  and  FFT  ; 

algorithms,  j 

The  lower  two  tables  are  both  included  to  show  that 

^ I 

there  is  no  discernible  difference  in  the  matrix  and  FFT  < I 


a 


T 
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approaches.  This  is  because  the  orthogonal  similarity  trans- 
formation on  the  real  symmetric  periodic  Whittaker  matrix, 
and  the  complex  unitary  equivalence  transformation  imple- 
mented via  the  FFT  must  multiply  to  the  same  matrix  except 
for  roundoff  error.  When  the  FFT  and  periodic  matrix 
algorithms  are  used  on  a perfectly  band-limited  function. 

Figure  8-11  shows  that  the  norm  of  roundoff  error,  ( |e|  |, 
from  the  matrix  technique  is  an  order  of  magnitude  greater 
than  for  the  FFT  algorithm. 

Figures  8-11  through  8-13  are  tabular  summaries  of 
data  from  Figures  8-1  through  8-10.  Figure  8-11  is  a 
summary  of  the  norms  of  the  original  data,  the  various  inter- 
polants  and  the  errors,  plus  a listing  of  the  program 
compilation  times  and  run  times.  All  the  data,  for  instance 
in  Figures  8-la,  b,  and  c,  were  produced  during  one  computer 
run  at  a cost  of  16.8  seconds  compile  time  and  11.19  seconds 
run  time  on  the  UNM  IBM  360-67  using  the  standard  Fortran 
IV-G  compiler.  Figure  8-12  has  scale  factors  applied 
(listed  in  heading  information  in  Figures  8-lc  to  8-lOc)  to 
the  maximum  and  minimum  errors  so  absolute  error  comparisons 
can  be  made.  Finally,  Figure  8-13  is  a comparison  of  the 
relative  maximum  errors  in  the  interpolants  and  is  the  most 
important  of  the  three  summaries. 

Figure  8-13  also  lists  a percentage  factor  for  how 
much  transient  signal  lies  outside  the  region  of  interpolation; 
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that  is,  given  that  each  data  vector  x behaves  inside  the 
region  of  interpolation  as  depicted  by  the  solid  curve  in 
Figures  8-la  to  8-lOa,  how  much  more  signal  must  there  be 
in  order  that  x come  from  a band-limited  process?  This 
amount  is  expressed  in  Figure  8-13  as  a percentage  of  | |x| | 
and  is  the  "a"  factor  derived  in  equation  7-31  times  100. 

Figure  8-13  summarizes  several  interesting  results. 
First,  and  obvious,  small  relative  error  implies  small  a. 

In  other  words,  if  we  properly  sample  a reasonably  band- 
limited  phenomena  and  use  2TW  seimples  in  the  transient 
interpolation  scheme,  then  we  should  expect  small  errors 
in  the  interpolants  and  little  significant  data  outside 
the  interpolation  interval.  This  result  is  clearly  demon- 
strated in  Signals  1,  2,  3,  4,  5,  8,  and  9 where  relative 
errors  are  less  than  .091  and  energy  outside  the  interpolation 
interval  (remember  spread  over  infinity)  is  less  than  43% 
of  the  energy  inside  the  interval.  Second,  FFT  (or  the 
equivalent  symmetric  periodic  matrix)  interpolation  reduces 
the  relative  error  when:  (a)  the  periodic  extension  of  the 
transient  phenomena  forms  a perfectly  band-limited  process; 
e.g.,  the  extension  of  the  single  pulse  of  a sinusoidal 

wave  form  in  Signal  1;  (b)  the  periodic  extension  of  the 

transient  phenomena  reduces  the  severity  of  a discontinuity; 
i.e.,  the  extension  of  the  half  cycle  sine  wave  plus  d.c. 

offset  in  Signal  3.  The  Figure  also  shows  that  FFT  inter- 
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polation  can  increase  the  error  when  the  phenomena  begins 
and  ends  at  significantly  different  levels  as  in  Signal  4. 
This  obtains  from  forcing  the  periodic  extension  and  thereby 
forcing  the  last  interpolating  point  to  be  midway  between 
the  discontinuity.  Finally,  and  surprisingly  so,  the  inter- 
polants  computed  by  the  selected  eigenvector  algorithm  are 
generally  better  in  the  sense  of  smaller  relative  errors  than 
those  using  the  full  transient  interpolation.  This  can  be 
seen  for  Signals  3,  4,  6,  8,  and  9.  A visual  improvement  is 
also  noticed  for  Signal  7 (Figure  8-7b)  where  the  oscillatory 
overshoot  of  full  transient  interpolation  is  significantly 
reduced  by  the  eigenf liter  techniques.  The  maximum  relative 
error,  though,  for  this  example  is  a little  larger  than  for 
full  transient  interpolation.  Actually,  to  this  author, 
every  example  used  visually  appears  - overall  - as  good  as 
or  better  when  using  the  selected  eigenvector  approach. 

This  assertation  is  verified  by  Figure  8-11  wherein  the  norms 
of  the  total  errors  using  selected  eigenvectors  are  less 
than  for  the  transient  errors  for  all  signals  except  1,  2, 
cind  5 and  approximately  equal  in  these  three  examples. 

8 . 3 Conclusions 

Our  main  intent  for  this  dissertation  was  to  describe 
interpolation  as  a matrix  process  and  to  present  algorithms 
which  implement  the  matrix  techniques  as  an  alternative  to 
simple  linear  ..  terpolation  algorithms.  Our  purpose  was  to 
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gain  additional  insight  into  what  interpolation  really  does 
to  the  engineering  band-limited  function. 

In  Chapters  IV,  V,  and  VI  we  did  describe  interpolation 
in  terms  of  the  Whittaker  matrix  processes.  We  were  able  to 
derive  closed  form  expressions  for  elements  of  the  periodic 
symmetric  matrix  in  terms  of  simple  cosecant  and  cotangent 
functions.  In  Chapter  V we  were  able  to  show  the  complete 
agreement  of  a symmetric  matrix  decomposition  description  of 
interpolation  and  the  Fast  Fourier  Transform  (FFT)  imple- 
mentation so  in  vogue  with  the  engineers.  In  Chapter  VI 
we  also  described  a new  matrix  algorithm  for  curve  fitting 
via  the  Eigen-Filter,  Cross-Correlation  Algorithm.  Here, 
we  showed  that  the  correlation  coefficient  is  the  indicator 
of  goodness  of  fit  of  a curve  to  data,  much  akin  to  the  way 
we  say  a minimum  sum  of  squared  residuals  is  a measure  of 
goodness  of  fit  of  a curve  in  the  least  squares  sense. 

The  new  insights  into  interpolation  are  implicit  in 
the  linear  algebra  interpretations  of  matrix  transformations 
on  linear  vector  spaces.  Given  2TW  samples  and  the  assur- 
ances of  the  data  gatherer  that  the  samples  are  not  aliased, 
we  can  now  view  interpolation  as  a mapping  of  the  data 
vector  of  2TW  components  into  an  interpolant  vector  of  2TW 
components.  Of  particular  importance  is  that  by  modifying 
the  transient  Whittaker  matrix  as  described  in  Chapter  VI, 
we  can  generate  an  orthogonal  similarity  transformation 
which  has  distinct  eigenvalues  and  a complete  set  of  unique 
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eigenvectors.  Any  function  in  the  space  of  2TW  dimension 
is  then  representable  as  a linear  combination  of  the  2TW 
basis  vectors  - the  eigenvectors.  We  used  10  rather  "nasty” 
functions,  in  the  Fourier  sense,  to  demonstrate  the  viability 
of  using  selected  eigenvectors  for  interpolation.  VJe  also 
showed  in  Chapter  VII  that  error  bounds  for  matrix  equations 
can  be  easily  and  clearly  established  in  terms  of  matrix  and 
vector  condition  nurabers  and  norms. 

As  an  alternative  to  linear  approximation  approaches  to 
interpolation,  we  presented  a whole  appendix  of  computer 
programs  to  implement  several  exact  interpolation  schemes. 

We  do  not  claim  our  code  is  optimum,  but  merely  that  it  does 
work  as  shown  by  the  numerous  examples  in  section  2 of  this 
chapter . 

In  conclusion,  we  point  out  that  any  finite  impulse 
response,  nonrecursive  digital  filter  can  be  formulated 
either  as  a transient  or  periodic  symmetric  matrix  process 
as  we  have  done.  Possibilities  for  future  work  abound:  by 
replacing  our  subroutine  MATRIX  with  one  that  generates  the 
appropriate  matrix  (or  its  inverse)  for  the  problem  at  hand, 
we  can  describe  filtering  or  deconvolution  (matrix  inverse) 
in  terms  of  similarity  transformations  and  matrix  norms. 
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PLOT  OF  PATRIX  INTERPOLATION 
HULTlPLIERa  0«tOOE  01  NlflNOa  0 
eigenvectors  used  in  THE  RECONSTRUCTION 
12345670 

9 10  11  12  13  14  IS  16 
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COLUMN  2 IS  THE  ERROR 
COLUMNS  PLOTTED#  1*  2t 


1 
2 

3 

4 

5 

6 

7 

8 
9 

10 
1 1 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 
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8 -la  Transient  Interpolation  and 
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PLOT  OP  MATRIX  INTERPOLATION 
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8-3b  Eigenfilter  and  FFT  Interpolation  - Signal  3 
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-4c  Tabular  Comparison  of  Interpolation 
Schemes  - Signal  4 


PLOT  OF  MATRIX  INTERPOLATION 
MULTIPLIER*  0«960E  90  NMINO*  0 
EIGENVECTORS  USED  IN  THE  RECONSTRUCTION 
12245678 

9 10  11  12  13  14  15  16 

COLUMN  1 IS  INTERPCLAT60  CATA 
COLUMN  2 IS  THE  ERRCR 
COLUMNS  PLOTTEOt  1*  2* 


LEFT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 

NUMBER  CF  PEAK-TO-PEAK  PRINT  CELLS*  27 


PLOT  CF  CORRELATICK  COEFFICIENTS 
MULTIPLIER*  0«100E  01  NMINO*  0 
COLUMNS  PLOTTED*  1* 


LEFT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 

NUMBER  OF  PEAK-TO-PEAK  PRINT  CELLS*  19 


8-5a  Transient  Interpolation  and 

Correlation  Coefficients  - Signal  5 
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PLOT  OF  MATRIX  INTERPOLATION 
NLLTlPLIERs  0«9R2E  00  N«InO«  0 
EIGENVECTORS  USED  IN  THE  RECONSTRUCTION 
00000000 

0011  00000 

COLUMN  I IS  INTERPCLATEO  CATA 
COLUMN  2 IS  THE  ERROR 
COLUMNS  PLOTTED*  1*  2* 


NUMBER  CF  PEAK-TO>PEAK  PRINT  CELLS*  27 


PLOT  OF  FFT  1 NTERPCLAT 1 CN 
MULTIPLIER*  0*960E  00  N«IN0«9S9 
COLUMN  1 IS  INTERPOL ATEO  CATA 
COLUMN  2 IS  THE  ERROR 
COLUMNS  PLOTTED*  1*  2* 


RIGHT  IS  POSITIVE 


ji 


NUMBER  CE  PEAK-TO-PEAK  PRINT  CELLS*  27 


8-5b  Eigenfilter  and  FFT  Interpolation  - Signal  5 
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PRINT  or  MATRIX  INTERPOLATION 
MULTIPLIER*  0.960e  00  N«lNO*  0 
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-5c  Tabular  Comparison  of  Interpolation 
Schemes  - Signal  5 


PLOT  OF  WATRI)<  INTCPPOLAT  tON 
PULTlPLlgPs  C*i30e  01  N«1N0«  0 

EIGENVCCTOPS  USED  IN  THE  PECO ^ STRUCT  I C K 
I 2 3 4 S 6 7 6 

^ 10  11  12  13  14  15  16 

COLUMN  1 IS  INTERPOLATED  DATA 
COLUMN  2 IS  the  error 
COLUMNS  plotted*  1«  2* 


LEFT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 


NUMBER  OF  PEAK-TO-PEAK 


PRINT  CELLS*  27 


PLOT  OF  CORRELATION  COEFFICIENTS 
multiplier*  C*103E  01  NMIND*  0 
COLUMNS  PLOTTED*  1* 


LEFT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 


NUMBER  OF  PEAK-*T0*PEAK  PRINT  CELLS*  10 


8-6a 


Transient  Interpolation  and 
Correlation  Coefficients  - Signal  6 
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PLOT  OF  PATRIX  INTEPPOLAT ICN 
MLLTlPLlCPs  0«126E  Ot  N«INO«  0 
EIGENVECTORS  USEO  IN  THE  RECONSTRUCT! CN 
OOOOSGOO 

C 0 11  12  13  0 15  0 

COLUMN  1 IS  INTERPCLATEO  CATA 
COLUMN  2 IS  THE  ERPCR 
COLUMNS  PLOTTED*  I*  2* 


LEFT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 


NUMBER  OF  PEAK*TO-PEAK  PRINT  CELLS*  27 


PLOT  OF  FFT  I NTERPCLAT I CN 
MULTIPLIER*  0«128E  01  N«IN0«9S9 
COLUMN  1 IS  INTERPCLATEO  CATA 
COLUMN  2 IS  THE  ERROR 
COLUMNS  PLOTTED*  1*  2* 


LEFT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 


NUMBER  CF  P6AK*ro*PEAK  PRINT  CELLS*  27 


8-6b  Eigenfilter  and  FFT  Interpolation  - Signal  6 


PKINT  CP  MATRIX  interpolation 
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-6c  Tabular  Comparison  of  Interpolation 
Schemes  - Signal  6 


PLOT  CP  MATRIX  interpolation 
MULTIPLIERS  0«117C  01  NRiNOs  0 
£ZC£NVCC70R$  USEO  IN  THE  RECOASTRUCT I ON 
00000070 

0 10  0000  00 

COLUMN  t IS  INTERPCLATEO  CATA 
COLUMN  2 IS  THE  ERROR 
COLUMNS  RLOTTEO*  1*  2* 
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LEPT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 


NUMBER  OF  PEAK-TO-PEAK  PRINT  CELLSa  27 


PLOT  OF  FFT  INTERPCLAT ICN 
MULTIPLIERS  0«128E  01  N«IN0a999 
COLUMN  1 IS  INTERPCLATEO  DATA 
COLUMN  2 IS  THE  ERROR 
COLUMNS  PLOTTED*  1*  2* 


LEFT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 

NUMBER  CF  PEAK«TO*PEAK  PRINT  CELLS*  27 


8-7b  Eigenfilter  and  FFT  Interpolation  - Signal  7 
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-7c  Tabular  Comparison  of  Interpolation 
Schemes  - Signal  7 


y 


PLOT  OF  MATRIX  INTERPOLATION 
multipliers  0*S266  00  NMINQs  0 
EIGENVECTORS  USED  IN  THE  RECO NSTRUCT I C K 
12  2 4 5 6 7 8 

9 10  II  12  13  14  IS  16 

COLUMN  1 IS  INTERPOLATED  DATA 
COLUMN  2 IS  THE  ERROR 
COLUMNS  PLOTTED*  1*  2* 


PLOT  OF  CORRELATICK  COEFFICIENTS 
MULTIPLIERS  0«100fi  01  NMlNOs  0 
COLUMNS  PLOTTED*  1* 
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LEFT  OF  CENTER  LINE  IS  NEGATIVE 
RIGHT  IS  POSITIVE 

NUMBER  CF  PEAK-rO-PEAK  PRINT  CELLS*  19 


8-8a  Transient  Interpolation  and 

Correlation  Coefficients  - Signal  8 
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-8c  Tabular  Comparison  of  Interpolation 
Schemes  - Signal  8 


PLOT  OF  FATPIX  INT EPPOLAT ION 
MLLTlPLIEPa  0*737E  00  N«1N0*  0 

FIG6NVEC70PS  USED  IN  THS  PECONSTPUCT f CK 
12  3 4 5 6 7 0 

9 1C  I 1 12  13  lA  IS  16 

COLUMN  I IS  INTERPOLATEO  DATA 
COLUMN  2 IS  THE  EPPOP 
COLUMNS  PLOTTED*  1*  2* 
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Figure  8-12.  Comparison  of  Absolute  Errors 
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CHAPTER  IX 


RECOMMENDATIONS  FOR  FUTURE  WORK 

We  set  out  in  this  work  to  propose  alternatives  to 
linear  interpolation  - these  we  have  demonstrated  in  previous 
chapters  through  the  Whittaker  matrix  processes.  We  end  our 
effort  by  proposing  several  directions  in  which  matrix  inter- 
polating techniques  might  be  extended.  As  is  most  often  the 
case  with  a little  knowledge  gained,  we  end  by  asking  more 
questions  than  we  answered;  however,  we  believe  the  problem 
areas  are  important  and  have  considered  each  in  some  detail. 

There  are  four  problems  that  we  would  like  to  pursue: 

(1)  extrapolation  outside  the  original  data  interval  using 
matrix  techniques;  (2)  inverse  interpolation  for  the  case 
when  the  Whittaker  matrix  is  singular;  (3)  calculating  the 
derivative  of  a seimple  set  at  each  sample  point  using  matrix 
techniques;  (4)  recursive  or  binary  interpolation  in  which 
the  matrix  equations  are  simplified.  We  formulate  and  discuss 
these  problems  in  the  following  sections. 

9 . 1 Extrapolation  Using  the  Transient  tVhittaker  Matrix 

Referring  back  to  Figure  4-1,  we  note  that  when  r = 0, 
f^  and  f^  for  the  two  choices  of  N are  really  "extrapolated" 
points.  This  is  inferred  in  the  sense  that  these  points  are 
outside  the  "original"  data.  We  pose  the  question  of  what 
happens  when  the  interpolants  and  the  extrapolant  themselves 
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are  interpolated.  When  N = 4 in  Figure  4-1,  we  would  expect 
this  second  interpolation  to  return  the  original  data 
and  plus  a new  extrapolated  point  (where  Xj^  is  located 
in  the  Figure) . Ignoring  errors  for  the  moment,  we  can 
continue  this  process  using  equation  4-8  as  follows: 

Define 


Then , 


f<°>  = X 


f(l)  ^ pf(0) 


f = P^f 


(9-1) 


(9-2) 


VThere  f is  a vector  and  the  bracketed  superscript  (i)  infers 
that  the  N element  of  the  vector  is  the  (i^“)  extrapolant. 

Equation  9-2  can  be  rewritten  in  terms  of  the  symmetric 
Whittaker  matrix 


or, 


g(N-i+l)  ^ ^(i) 


(9-3) 


(9-4) 


Using  the  orthogonal  similarity  transformation  of  Chapter  VI 
^ (QrQ'^)^g^°^  (9-5) 


or 


(9-6) 
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Thus,  it  seems,  extrapolation  could  be  implemented  by  the 
programs  in  the  Appendix  upon  raising  the  eigenvalues  of  S 
to  the  appropriate  power.  We  now  show,  however,  that  this 
is  unwise. 

Using  the  approach  in  Chapter  VII,  we  try  to  find  the 
original  data  vector  from  our  N-times  interpolated  inter- 
polants;  i.e.. 


* 

X 


„-N  (n) 
S g 


(9-7) 


This  implies  that  after  N successive  interpolations  using 
equation  9-6  we  cannot  recover  x but,  rather,  some  "nearby" 
vector  X*.  Then,  proceeding  as  in  Section  5 of  Chapter  VII 
we  find 


= S^(x  - X*)  = S^r 


(9-8) 


where  e^^^  is  the  error  vector  for  g^^^  and  r is  the  residual 
vector.  Then 


(N) 


,Ni 


,-N, 


(N) 


= [k(S)] 


N r 


(9-9) 


or, 


(N) 


“(nT 


y N 
max  I 


'min 


1^1 


(9-10) 


We  note  that  N = 16,  produces  a magnification  factor  on  the 


order  of  10 Thus  it  seems  that  equation  9-6  is  unsatis- 


factory.  However,  if  we  expand  g 


(1) 

t 


I 


and  since  we  know  the  x^,  we  can  determine  a 


(9-11) 


(9-12) 


where  Ag^  must  be  estimated,  perhaps  from  equation  7-31. 
We  can  then  form  a corrected 


g(l)  = g(l)  - 


(9-13) 


and  continue  as 


g(2)  = QrVg<^^ 


(9-14) 


(i) 


at  each  step  in  terms  of  the  previous 


By  computing  the  g 
^ ( i — 1 ) 

g , It  may  be  possible  to  significantly  reduce  the 
(N) 

errors  in  g . Needed,  of  course,  are  ways  to  estimate 
Ag^^  and  a rigorous  error  analysis  of  the  performance  of 
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equation  9-6  using  corrected  data. 

9 . 2 Inverse  Interpolation 

The  inverse  problem  is  simply  stated  as:  given  a 
vector  g of  interpolants  (perhaps  computed  by  some  linear 
approximation  scheme) , how  can  we  find  the  original  data 
vector  X?  If  we  consider  the  interpolants  as  the  discrete 
output  from  a discrete  filter,  we  ask  how  do  we  find  the 
discrete  input  samples.  From  our  discussions  in  Chapter  VII 
concerning  errors,  we  know  in  general  that  computing 
X = S ^g  returns  not  x,  but  rather,  x*,  some  "nearby"  vector. 
Suppose  though,  that  S ^ does  not  exist.  This  is  exactly 
the  problem  in  the  periodic  Whittaker  matrix  with  even 
ordered  dimensions.  We  might  ask,  then,  what  this  means  in 
terms  of  our  error  bounds,  especially  since  g = Sx  is  a 
perfectly  valid  way  of  finding  g given  x. 

Several  important  results  from  Linear  Algebra  help  put 
this  problem  into  perspective  [10,  p.  49].  First,  Sx  = g 
has  NO  solution  x for  certain  g.  Also,  there  are  non-zero 
vectors  x which  satisfy  Sx  = 0 . Finally,  there  are  non- 
unique solution  vectors  x to  Sx  = g. 

The  first  of  these  results  is  the  most  interesting 
because  it  implies  that  we  can  produce  interpolants  from 
some  processes  which  cannot  possibly  come  from  interpolating 
with  the  Whittaker  theory.  In  fact,  [10,  p.  48]  any 
vector  g which  is  not  a linear  combination  of  the  columns 
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of  S cannot  be  produced  by  Whittaker  interpolation.  This 
result  obtains  from  the  well  known  theorem  that  when  the 
vector  g is  appended  to  the  columns  of  S,  the  rank  of  the 
resulting  S must  not  change.  If  g is  not  a linear  combination 
of  the  columns  of  S,  then  rank  S will  change  and  thus  there 
is  no  solution  of  Sx  = g.  This  all  implies  that  even  ordered 
Whittaker  interpolation  excludes  certain  sample  sets. 

As  to  the  last  two  results,  we  know  that  any  vector 


is  a solution  of  Sx  = 0 when  order  of  S is  even.  This  is 
obvious  from  the  form  of  equations  5-14  and  5-16  wherein  the 
sum  of  the  column  vectors  of  S is  zero  if  the  sign  of  each 
column  is  varied  by  equation  9-15.  Thus,  for  any  solution 
X of  Sx  = g,  we  know  x + x is  also  a solution.  The  vector  x 
is  simply  the  sample  set  from  any  periodic  function  whose 
frequency  is  exactly  equal  to  the  Nyquist  rate. 

Further  study  is  needed  of  the  error  bounds  for  even 
ordered  interpolation.  The  mechanics  of  bounding  errors  in 
g = Sx  in  terms  of  the  actual  equation  implementing  the 
matrix  products  is  well  established  in  numerical  analysis. 
Also,  the  numerical  errors  in  computing  x (e.g.,  by  Gaussian 


elimination)  are  also  established.  These  results  can  be 
adapted  to  the  special  forms  of  the  Whittaker  matrix  equations. 


9 . 3 Computing  the  Derivative 

One  of  many  possible  numerical  formulas  for  estimating 
the  derivative  of  a function  using  its  sample  set  is  [28, 
p.  98] 


X . , . - X . , 

• j.  i+h  i-h 

i 2h 


(9-16) 


where  h = j the  sampling  interval  if  samples  from  the  mid- 
point Whittaker  process  are  to  be  used.  Our  problem  here  is 
to  outline  a study  of  calculating  the  derivative  at  all  the 
Scimple  points  using  some  form  of  the  above  equation. 

First,  recognize  that  the  periodic  interpolating 
equation 


f = Px 


(9-17) 


produces  a vector  of  interpolants  at  the  half  steps.  Further- 
more, if  we  should  multiply  f by  a permutation  matrix  of  the 
form 


(9-18) 
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we  shift  all  elements  of  f to  i - j.  Then  we  can  write 


. . f - Tf  (P  - TP) 

^ “ 2h  2h  ^ 


(9-19) 


which  is  the  matrix  equation  for  computing  the  derivative 
at  all  the  original  data  points  x^ . 

From  Chapter  V,  P is  cyclic;  then 


^ ^ (VAV*  - TVAV*) 
2h  ^ 


(9-20) 


But  T is  cyclic  too.  From  Grey  [13,  pp.  16-21] , all  cyclic 
matrices  have  the  same  eigenvectors;  therefore 


. . (VAV*  - vyv*vAv*) 

^ 2h  ^ 


V[A(I  - 7)]V* 
2h 


(9-21) 


(9-22) 


As  shown  in  Chapter  V,  equation  9-22  is  implementable  via 
the  FFT  algorithm.  We  also  know  that  the  algorithm's 
numerical  problems  are  expressable  as 


_1_  I \nax  I 
" 2h  'a 


max 


(9-23) 


min 


where  ^ is  the  error  in  the  derivatives,  the  are  eigen- 
values of  (I  - V),  and  r is  x*  - x with  x*  computed  from 
the  inverse  of  equation  9-19.  We  have  already  shown  that 


W, 


min 


= 1 for  the  periodic  casj  (N  * 9) , so  the 
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stability  of  the  derivative  algorithm  is 


(9-24) 


We  see  that  small  h and/or  eigenvalues  of  T close  to  1 

(1  - tj^  -*•  0)  cause  the  error  bound  to  blow  up.  In  fact, 

for  the  particular  form  of  T chosen  (equation  9-18)  we 

know  that  the  t.  are  the  N roots  of  1.  Therefore,  a . = 0 

1 min 

and  equation  9-24  is  unbounded. 

We  propose  that  other  derivative  algcilthms  be  in- 
vestigated to  find  those  with  finite  (well  conditioned) 
error  bounds.  Also,  the  derivative  formulae  should  be 
extended  to  the  transient  Whittaker  matrix. 


9 . 4 Recursive  (Binary)  Interpolation 

Again,  let  the  superscript  represent  the  number  of  data 
points  used  in  an  interpolation  scheme.  Since  S is  always 
square,  N also  represents  the  number  of  interpolants  pro- 
duced. Then,  we  can  form  the  following  sequences 


N N 

g = S X 


(9-25) 


2N  N , N 

X = T^g  + T2X 


(9-26) 


N 


where  Tj^  and  T2  are  permutation  matrices  which  allow  g and 
N 

X to  be  added  together;  for  example,  when  N = 2 
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0 0 


0 1 


0 0 


1 0 


1 0 


0 0 


0 1 


0 0 


- g- 


(9-27) 


Continuing, 


2N  „2N  2N 
g = S X 


^2N,  N N ^ „N  N. 
S (Tj^g  + T2X  ) 


s2N(tJs^  + T^)x^ 


4N  „2N  2N  „2N  2N 
X = g + T2  X 


g4N  ^ 34N^4N  ^ 34N(  2Ng2N^^2N^2Nj 


(9-28) 

(9-29) 


= + T2^)  (T^S^  + T2)x^  (9-30) 


ginN  ^ gmN^mN  ^ ^mN  ^ + Tf)]x^  (9-31) 


where  m=  2,  4,  8,  16,  ...  The  dimensions  of  the  resulting 


matrix  equation  can  be  shown  as 


X 


(9-32) 


mNxl 


mN  X mN 


mNxN 
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All  the  problems  attacked  in  this  dissertation  now  have 
analogues  in  the  non-square  system  expressed  by  equation 
9-32.  One  potentially  fruitful  approach  to  their  solution 
might  be  to  investigate  the  behavior  of  the  "singular  values" 
[2],  [33]  of  the  system  using  singular  value  decomposition 

(SVD) . This  approach  expresses 


g = 


(9-33) 


where  U and  V are  orthogonal  matrices  (different  dimensions) 
and  is  a non-square  matrix  of  the  form 


(9-34) 


mNxN 


T T 

with  the  the  square  roots  of  the  eigenvalues  of  (H  S sn) . 

T 

We  note  that  for  the  nonsingular  periodic  S,  N=9,  SS=I, 

and  that  the  are  then  the  square  roots  of  the  eigenvalues 
T 

of  (n  IT).  It  should  be  possible  to  formulate  a Singular- 
Filter,  Cross-Correlation  Algorithm  as  an  analogue  to  the 
Eigen-Filter,  Cross-Correlation  Algorithm  of  Chapter  VI. 
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